#---------------------------------#
# Libraries & Settings ####
#---------------------------------#


#Uncomment to install relevant packages
# install.packages("plyr")
# install.packages("zoo")
# install.packages("texreg")
# install.packages("readxl")
# install.packages("lfe")
# install.packages("dplyr")
# install.packages("R2HTML")
# install.packages("ggplot2")
# install.packages("Hmisc")
# install.packages("RColorBrewer")
# install.packages("devEMF")
# install.packages("gridExtra")


library("plyr",include.only = "rbind.fill")
library("zoo")
library("texreg")
library("readxl")
library("lfe")
library("dplyr")
library("R2HTML")
library("ggplot2")
library("Hmisc")
library(RColorBrewer)
library(devEMF)
library(gridExtra)


spectral  <-  brewer.pal(11, "Spectral")
darkcols  <-  c(brewer.pal(8, "Dark2"),brewer.pal(8, "Paired"))

setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
setwd("..")


#---------------------------------#
# Functions ####
#---------------------------------#
logodds100 <- function(x)
{
  logodds <- log(x/(100-x))
  return(logodds)
}
logodds <- function(x)
{
  logodds <- log(x/(1-x))
  return(logodds)
}
oddsratio <- function(x,y)
{
  or <- (x/(1-x))/(y/(1-y))
  return(or)
}
id <- function(x)
{
  id <- x*1
  return(id)
}
simplelag <-function(x,by=NULL,mylag=1,outside=NA)
{
  myend<-length(x)- mylag
  if (!is.null(by)) {
    lby<-c(replicate(mylag,""),as.character(by[1:myend]))
    y0<-c(replicate(mylag,outside),x[1:myend])
    y<-ifelse(as.character(by)==lby,y0,outside)
    
  }
  else {
    y<-c(replicate(mylag,outside),x[1:myend])
  }
}

mygraph_f  <- function(myvar,mytitle,my_y,mylog,ratio="logodds",legpos="topleft",average=3,myyaxt="s",
                       my_mar=c(2.5,4,2.5,1),my_inset=c(0,0.1),my_legend="y",myxaxt="s",my_ylim=NULL,my_xlim=NULL) {    
  
  m <- get(substr(deparse(substitute(myvar)),1,-1+regexpr("\\$",deparse(substitute(myvar)))))
  m$dv <- myvar
  
  g <- as.data.frame(tapply(myvar,list(m$year,m$country),mean))
  
  
  if (ratio=="OR" & myyaxt=="n") {myfunc <- logodds100} else  {myfunc <- id}
  
  g$empty <- NA
  delcol <- NULL
  for (i in 1:ncol(g)) {
    if(sum(is.na(g[,i]))>=nrow(g)-1) { delcol  <- c(delcol,i)}
  }
  g <- g[,-delcol]
  if (ratio=="OR") {
    g <- 100*g
    percent <- "%"
  } else {percent <- ""}
  g$year <- rownames(g)
  
  year <- data.frame(c(min(g$year):max(g$year)))
  colnames(year) <- "year"
  rownames(year) <- year$year
  g <- merge(year,g,by="year",all.x=TRUE)
  rownames(g) <- year$year
  g2 <- g
  for (i in 1:ncol(g)) {
    if(sum(is.na(g[,i]))<nrow(g)-1) { g2[,i]  <- na.approx(g[,i],rownames(g),na.rm=FALSE)}
  }
  g2dif <- g2[2:nrow(g2),]-g2[1:nrow(g2)-1,]
  
  if (average==3)
  { 
    h <- (g2+rbind(rep(NA,ncol(g2)),g2[1:nrow(g2)-1,])+rbind(g2[2:nrow(g2),],rep(NA,ncol(g2))))/3
  } else
  {
    h <- g2
  }
  h2 <- h
  h2dif <- h2[2:nrow(h2),]-h2[1:nrow(h2)-1,]
  
  y <- (1-is.na(h))*h$year
  y[y==0] <- NA
  y <- as.data.frame(y)
  mystart <- sapply(y,min,na.rm=TRUE)
  myend <- sapply(y,max,na.rm=TRUE)
  first0 <- (is.na(h)==FALSE & is.na(rbind(rep(NA,ncol(h)),h[1:nrow(h)-1,])==TRUE))*h
  first0[first0==0] <- NA
  first <- sapply(first0,max,na.rm=TRUE)
  last0 <- (is.na(h)==FALSE & is.na(rbind(h[2:nrow(h),],rep(NA,ncol(h))))==TRUE)*h
  last0[last0==0] <- NA
  last <- sapply(last0,max,na.rm=TRUE)
  
  
  if (ratio=="OR") 
  {
    my_ratio <- exp(log(oddsratio(last/100,first/100))/(myend-mystart))
  } else if (ratio=="logodds") 
  {
    
    my_ratio <- exp((last-first)/(myend-mystart))
    
  }  else if (ratio=="simple") 
  {
    
    my_ratio <- exp(log(last/first)/(myend-mystart))
    
  }
  
  yratio <- round(my_ratio,4)
  ylogratio <- round(log(my_ratio),4)
  gr_rate <- round((my_ratio-1)*100,2)
  
  
  g3 <- data.frame(year,rowMeans(h2[,-1],na.rm=TRUE))
  colnames(g3)[2] <- "mean_or"
  g3$ldif <- c(NA,rowMeans(h2dif[,-1],na.rm=TRUE))
  g3$nbval <- ncol(h2[,-1])-rowSums(is.na(h2[,-1]))
  g3$maxnbval <- max(g3$nbval)
  g3$atmax <- (g3$nbval>=g3$maxnbval)*1
  g3$aftermax <- (cumsum(g3$atmax)>=1)*1
  g3$meanco <- ifelse(g3$atmax,g3$mean_or,NA)
  
  for (i in 1:nrow(g3)) {
    if(is.na(g3$meanco[i])==TRUE & g3$aftermax[i]==1)  {g3$meanco[i] <- g3$meanco[i-1]+g3$ldif[i]}
  }
  g3 <- g3[order(rev(sort(g3$year))),]
  for (i in 1:nrow(g3)) {
    if(is.na(g3$meanco[i])==TRUE & g3$aftermax[i]==0)  {g3$meanco[i] <- g3$meanco[i-1]-g3$ldif[i-1]}
  }
  g3 <- g3[order(g3$year),]
  g3$meanco
  g3$meanco2 <- ifelse(g3$nbval>2,g3$meanco,NA)
  
  bottom <- min(g3$year[g3$nbval>2])
  top <- max(g3$year[g3$nbval>2])
  
  if (ratio=="OR") {
    r <- (g3$meanco[g3$year==top]/(100-g3$meanco[g3$year==top]))/(g3$meanco[g3$year==bottom]/(100-g3$meanco[g3$year==bottom]))
    r <- (exp(log(r)/(top-bottom))-1)*100
  }  else if (ratio=="logodds") {
    r <- g3$meanco[g3$year==top]-g3$meanco[g3$year==bottom]
    r <- (exp(r/(top-bottom))-1)*100
    if (myyaxt=="n") {g3$meanco <- exp(g3$meanco)}  
  } else if (ratio=="simple") {
    r <- g3$meanco[g3$year==top]/g3$meanco[g3$year==bottom]
    r <- (exp(log(r)/(top-bottom))-1)*100
  }
  
  mycol <- as.vector(na.omit(t(ctrycolor[1,which(names(ctrycolor) %in% colnames(g2[-1]))])))
  mylty <- as.vector(as.numeric(na.omit(t(ctrycolor[2,which(names(ctrycolor) %in% colnames(g2[-1]))]))))
  mypch <- as.vector(as.numeric(na.omit(t(ctrycolor[3,which(names(ctrycolor) %in% colnames(g2[-1]))]))))
  
  par(xpd=FALSE, mar=my_mar)
  # par(oma=c(0,0,0,0))
  matplot(g$year,myfunc(g[,2:ncol(g)]),bg=mycol,pch=mypch,col=mycol,ylab=my_y,xlab="",cex=0.7,las = 2, cex.lab=0.7, 
          cex.axis=0.7, log=mylog,yaxt=myyaxt,xaxt=myxaxt,xlim=my_xlim,
          ylim=my_ylim)
  if (myxaxt!="n") 
  {  
    axis(1, at = c(1971:1974,1976:1979,1981:1984,1986:1989,1991:1994,1996:1999,2001:2004,2006:2009,2011:2014,2016:2019,6),las=2,cex.axis=0.7)
  }
  matlines(g2$year,myfunc(cbind(g2[,2:ncol(g)])),lty=mylty,col=mycol)
  matlines(g3$year,myfunc(g3$meanco2),lwd = 3,col="black")
  if (ratio=="logodds" & myyaxt=="n")   
  {
    scale <- c(1/10,1/9,1/8,1/7,1/6,1/5,1/4,1/3,1/2.5,1/2,1/1.8,1/1.6,1/1.4,1/1.2,1,1.2,1.4,1.6,1.8,2,2.5,3,4,5,6,7,8,9,10)
    scalelabel <- c("/10","/9","/8","/7","/6","/5","/4","/3","/2.5","/2",
                    "/1.8","/1.6","/1.4","/1.2","1","× 1.2","× 1.4","× 1.6","× 1.8",
                    "× 2","× 2.5","× 3","× 4","× 5","× 6","× 7","× 8","× 9","× 10")
    
    axis(2,log(scale),label =scalelabel, las=2,cex.axis=0.8)  
  } else if (ratio=="OR" & myyaxt=="n")   
  {
    prop <- c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,
              1,2,3,4,5,6,7,8,9,
              10,12,14,16,18,20,25,30,35,40,45,50,55,60,65,70,75,80,82,84,86,88,90,
              91,92,93,94,95,96,97,98,99,
              99.1,99.2,99.3,99.4,99.5,99.6,99.7,99.8,99.9)
    axis(2,logodds100(prop),label = paste(prop,"%"), las=2,cex.axis=0.8)
  } 
  
  par(xpd=TRUE)
  if (my_legend!="n")
  {    
    legend(legpos, 
           inset=my_inset,
           legend=c(colnames(g[,-1]),"Adj. Mean"),pch =c(mypch,NA), 
           pt.bg=c(mycol,"black"),bg=NA,col=c(mycol,"black"),cex=0.8,ncol=7,lty=c(mylty,1),lwd = c(rep(1,ncol(g)-1),3),
           text.col = c(mycol,"black"),box.col="gray",box.lty=0)
  }
  c(g3$meanco[g3$year==bottom], g3$meanco[g3$year==top], g3$meanco[g3$year==top]/g3$meanco[g3$year==bottom])
  if (mytitle!="") {
    title(main=mytitle,sub=paste(ncol(g)-1,"countries.",bottom,":",round(g3$meanco[g3$year==bottom],1),percent,";",top,":", round(g3$meanco[g3$year==top],1),percent,";",
                                 "Yearly growth rate :",round(r,1),"%"))
  }
  if (r>0){plus <- "+"
  }  else {plus <- ""}
  mydelta <- paste("\n",plus,round(r,1),"% / year")
  
  if (ratio=="logodds") {  
    percent2 <- ""
    myfunc2 <- exp
    cross <- "×"
  } else {
    percent2 <- "%"
    myfunc2 <- id
    cross <- ""
  }
  
  text(g3$year[g3$year %in% c(bottom,top)],myfunc(g3$meanco2[g3$year %in% c(bottom,top)]),
       labels=paste(g3$year[g3$year %in% c(bottom,top)],"\n",
                    cross,round(myfunc2(g3$meanco2[g3$year %in% c(bottom,top)]),1),percent2,sep=""),
       pos=c(2,4),cex= 0.7)
  
  des <- data.frame(mystart,first,myend,last,yratio,ylogratio,gr_rate)
  des$country <- rownames(des)
  des1 <- des[1,]
  des1$country <- "Adj. Mean"
  des1$mystart <- bottom
  des1$first <- g3$meanco2[g3$year %in% c(bottom)]
  des1$myend <- top
  des1$last <- g3$meanco2[g3$year %in% c(top)]
  des1$gr_rate <- round(r,2)
  des1$yratio <- 1+r/100
  des1$ylogratio <- log(des1$yratio)
  des <- rbind.fill(des,des1)
  
  rownames(des) <- des$country
  
  
  
  print(des)
  
  if (ratio=="OR")
  {
    ltr <- lm(I(100*logodds(dv))~country+country:year,data=m[(m$country %in% "France" & m$year==1993)==F
                                                             & (m$country %in% "Denmark" & m$year==1994)==F,])
    ltr_all <- lm(I(100*logodds(dv))~country+year,data=m[(m$country %in% "France" & m$year==1993)==F
                                                         & (m$country %in% "Denmark" & m$year==1994)==F,])
    ltr_all2 <- lm(I(100*logodds100(meanco2))~year,data=g3)
  } else if (ratio=="simple")
  {
    ltr <- lm(I(100*log(dv))~country+country:year,data=m[(m$country %in% "France" & m$year==1993)==F
                                                         & (m$country %in% "Denmark" & m$year==1994)==F,])
    ltr_all <- lm(I(100*log(dv))~country+year,data=m[(m$country %in% "France" & m$year==1993)==F
                                                     & (m$country %in% "Denmark" & m$year==1994)==F,])
    ltr_all2 <- lm(I(100*log(meanco2))~year,data=g3)
  } else if (ratio=="logodds")
  {
    ltr <- lm(I(100*dv)~country+country:year,data=m[(m$country %in% "France" & m$year==1993)==F
                                                    & (m$country %in% "Denmark" & m$year==1994)==F,])
    ltr_all <- lm(I(100*dv)~country+year,data=m[(m$country %in% "France" & m$year==1993)==F
                                                & (m$country %in% "Denmark" & m$year==1994)==F,])
    ltr_all2 <- lm(I(100*meanco2)~year,data=g3)
  }
  print(screenreg(list(ltr,ltr_all,ltr_all2),digits=4))
  
  return(list(m=m,g=g,g2=g2,g3=g3,h2=h2,mycol=mycol,des=des
              ,ltr=ltr,ltr_all=ltr_all,ltr_all2=ltr_all2,mydelta=mydelta
  ))
  
}


mygraph2  <- function(myvar,mytitle,my_y,mylog,ratio="logodds",legpos="topleft",average=3,myyaxt="s") {    
  
  m <- get(substr(deparse(substitute(myvar)),1,-1+regexpr("\\$",deparse(substitute(myvar)))))
  # myvar2 <- get(substr(deparse(substitute(myvar)),1+regexpr("\\$",deparse(substitute(m))),deparse(substitute(m))))
  # m$dv2 <- m$myvar2
  m$dv <- myvar
  
  g <- as.data.frame(tapply(myvar,list(m$year,m$country),mean))
  
  
  #supprimer colonnes avec donnees vides
  if (ratio=="OR" & myyaxt=="n") {myfunc <- logodds100} else  {myfunc <- id}
  
  g$empty <- NA
  delcol <- NULL
  for (i in 1:ncol(g)) {
    if(sum(is.na(g[,i]))>=nrow(g)-1) { delcol  <- c(delcol,i)}
  }
  g <- g[,-delcol]
  if (ratio=="OR") {
    g <- 100*g
    percent <- "%"
  } else {percent <- ""}
  g$year <- rownames(g)
  
  year <- data.frame(c(min(g$year):max(g$year)))
  colnames(year) <- "year"
  rownames(year) <- year$year
  g <- merge(year,g,by="year",all.x=TRUE)
  rownames(g) <- year$year
  g2 <- g
  for (i in 1:ncol(g)) {
    if(sum(is.na(g[,i]))<nrow(g)-1) { g2[,i]  <- na.approx(g[,i],rownames(g),na.rm=FALSE)}
  }
  g2dif <- g2[2:nrow(g2),]-g2[1:nrow(g2)-1,]
  
  if (average==3)
  { 
    h <- (g2+rbind(rep(NA,ncol(g2)),g2[1:nrow(g2)-1,])+rbind(g2[2:nrow(g2),],rep(NA,ncol(g2))))/3
  } else
  {
    h <- g2
  }
  h2 <- h
  h2dif <- h2[2:nrow(h2),]-h2[1:nrow(h2)-1,]
  
  y <- (1-is.na(h))*h$year
  y[y==0] <- NA
  y <- as.data.frame(y)
  mystart <- sapply(y,min,na.rm=TRUE)
  myend <- sapply(y,max,na.rm=TRUE)
  first0 <- (is.na(h)==FALSE & is.na(rbind(rep(NA,ncol(h)),h[1:nrow(h)-1,])==TRUE))*h
  first0[first0==0] <- NA
  first <- sapply(first0,max,na.rm=TRUE)
  last0 <- (is.na(h)==FALSE & is.na(rbind(h[2:nrow(h),],rep(NA,ncol(h))))==TRUE)*h
  last0[last0==0] <- NA
  last <- sapply(last0,max,na.rm=TRUE)
  
  
  
  if (ratio=="OR") 
  {
    my_ratio <- exp(log(oddsratio(last/100,first/100))/(myend-mystart))
  } else if (ratio=="logodds") 
  {
    
    my_ratio <- exp((last-first)/(myend-mystart))
    
  }  else if (ratio=="simple") 
  {
    
    my_ratio <- exp(log(last/first)/(myend-mystart))
    
  }
  
  yratio <- round(my_ratio,4)
  ylogratio <- round(log(my_ratio),4)
  gr_rate <- round((my_ratio-1)*100,2)
  
  g3 <- data.frame(year,rowMeans(h2[,-1],na.rm=TRUE))
  colnames(g3)[2] <- "mean_or"
  g3$ldif <- c(NA,rowMeans(h2dif[,-1],na.rm=TRUE))
  #g3$fdif <- c(rowMeans(g2dif[,-1],na.rm=TRUE),NA)
  g3$nbval <- ncol(h2[,-1])-rowSums(is.na(h2[,-1]))
  g3$maxnbval <- max(g3$nbval)
  g3$atmax <- (g3$nbval>=g3$maxnbval)*1
  g3$aftermax <- (cumsum(g3$atmax)>=1)*1
  g3$meanco <- ifelse(g3$atmax,g3$mean_or,NA)
  
  for (i in 1:nrow(g3)) {
    if(is.na(g3$meanco[i])==TRUE & g3$aftermax[i]==1)  {g3$meanco[i] <- g3$meanco[i-1]+g3$ldif[i]}
  }
  g3 <- g3[order(rev(sort(g3$year))),]
  for (i in 1:nrow(g3)) {
    if(is.na(g3$meanco[i])==TRUE & g3$aftermax[i]==0)  {g3$meanco[i] <- g3$meanco[i-1]-g3$ldif[i-1]}
  }
  g3 <- g3[order(g3$year),]
  g3$meanco
  g3$meanco2 <- ifelse(g3$nbval>2,g3$meanco,NA)
  
  bottom <- min(g3$year[g3$nbval>2])
  top <- max(g3$year[g3$nbval>2])
  
  if (ratio=="OR") {
    r <- (g3$meanco[g3$year==top]/(100-g3$meanco[g3$year==top]))/(g3$meanco[g3$year==bottom]/(100-g3$meanco[g3$year==bottom]))
    r <- (exp(log(r)/(top-bottom))-1)*100
  }  else if (ratio=="logodds") {
    r <- g3$meanco[g3$year==top]-g3$meanco[g3$year==bottom]
    r <- (exp(r/(top-bottom))-1)*100
    if (myyaxt=="n") {g3$meanco <- exp(g3$meanco)}  
  } else if (ratio=="simple") {
    r <- g3$meanco[g3$year==top]/g3$meanco[g3$year==bottom]
    r <- (exp(log(r)/(top-bottom))-1)*100
  }
  #mycol <- darkcols[1:ncol(g)-1]
  mycol <- as.vector(na.omit(t(ctrycolor[1,which(names(ctrycolor) %in% colnames(g2[-1]))])))
  mylty <- as.vector(as.numeric(na.omit(t(ctrycolor[2,which(names(ctrycolor) %in% colnames(g2[-1]))]))))
  mypch <- as.vector(as.numeric(na.omit(t(ctrycolor[3,which(names(ctrycolor) %in% colnames(g2[-1]))]))))
  
  # par(xpd=FALSE)
  # par(xpd = NA, mar = c(4, 4, 1, 1))
  # par(oma=c(0,0,0,0))
  matplot(g$year,myfunc(g[,2:ncol(g)]),pch=mypch,col=mycol,bg=mycol,ylab=my_y,xlab="",cex=0.7,las = 2, cex.lab=0.7, 
          cex.axis=0.7, log=mylog, yaxt=myyaxt)
  axis(1, at = c(191,68:191,69,1971:1974,1976:1979,1981:1984,1986:1989,1991:1994,1996:1999,2001:2004,2006:2009,2011:2014,2016:2019,6),las=2,cex.axis=0.7)
  matlines(g2$year,myfunc(cbind(g2[,2:ncol(g)])),lty=mylty,col=mycol)
  matlines(g3$year,myfunc(g3$meanco2),lwd = 3,col="black")
  if (ratio=="logodds" & myyaxt=="n")   
  {
    scale <- c(1/10,1/9,1/8,1/7,1/6,1/5,1/4,1/3,1/2.5,1/2,1/1.8,1/1.6,1/1.4,1/1.2,1,1.2,1.4,1.6,1.8,2,2.5,3,4,5,6,7,8,9,10)
    scalelabel <- c("/10","/9","/8","/7","/6","/5","/4","/3","/2.5","/2",
                    "/1.8","/1.6","/1.4","/1.2","1","× 1.2","× 1.4","× 1.6","× 1.8",
                    "× 2","× 2.5","× 3","× 4","× 5","× 6","× 7","× 8","× 9","× 10")
    axis(2,log(scale),label =scalelabel, las=2,cex.axis=0.8)  
  } else if (ratio=="OR" & myyaxt=="n")   
  {
    prop <- c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,
              1,2,3,4,5,6,7,8,9,
              10,12,14,16,18,20,25,30,35,40,45,50,55,60,65,70,75,80,82,84,86,88,90,
              91,92,93,94,95,96,97,98,99,
              99.1,99.2,99.3,99.4,99.5,99.6,99.7,99.8,99.9)
    axis(2,logodds100(prop),label = paste(prop,"%"), las=2,cex.axis=0.8)
  } 
  
  # par(xpd=TRUE)
  legend(legpos, 
         #inset=c(0,-0.20),
         #inset=.01, 
         legend=c(colnames(g[,-1]),"Adj. Mean"),pch =c(mypch,NA), 
         col=c(mycol,"black"),cex=0.7,ncol=3,lty=c(mylty,1),lwd = c(rep(1,ncol(g)-1),3),
         text.col = c(mycol,"black"),pt.bg=c(mycol,"black"),bg=NA,box.col="gray",box.lty=0)
  c(g3$meanco[g3$year==bottom], g3$meanco[g3$year==top], g3$meanco[g3$year==top]/g3$meanco[g3$year==bottom])
  title(main=mytitle, 
        sub=paste(ncol(g)-1,"countries.",bottom,":",round(g3$meanco[g3$year==bottom],1),percent,";",top,":", round(g3$meanco[g3$year==top],1),percent,";",
                  "Yearly growth rate :",round(r,1),"%"))
  
  
  des <- data.frame(mystart,first,myend,last,yratio,ylogratio,gr_rate)
  des$country <- rownames(des)
  des1 <- des[1,]
  des1$country <- "Adj. Mean"
  des1$mystart <- bottom
  des1$first <- g3$meanco2[g3$year %in% c(bottom)]
  des1$myend <- top
  des1$last <- g3$meanco2[g3$year %in% c(top)]
  des1$gr_rate <- round(r,2)
  des1$yratio <- 1+r/100
  des1$ylogratio <- log(des1$yratio)
  des <- rbind.fill(des,des1)
  
  rownames(des) <- des$country
  print(des)
  
  if (ratio=="OR")
  {
    ltr <- lm(I(100*logodds(dv))~country+country:year,data=m[(m$country %in% "France" & m$year==1993)==F
                                                             & (m$country %in% "Denmark" & m$year==1994)==F,])
    ltr_all <- lm(I(100*logodds(dv))~country+year,data=m[(m$country %in% "France" & m$year==1993)==F
                                                         & (m$country %in% "Denmark" & m$year==1994)==F,])
    ltr_all2 <- lm(I(100*logodds100(meanco2))~year,data=g3)
  } else if (ratio=="simple")
  {
    ltr <- lm(I(100*log(dv))~country+country:year,data=m[(m$country %in% "France" & m$year==1993)==F
                                                         & (m$country %in% "Denmark" & m$year==1994)==F,])
    ltr_all <- lm(I(100*log(dv))~country+year,data=m[(m$country %in% "France" & m$year==1993)==F
                                                     & (m$country %in% "Denmark" & m$year==1994)==F,])
    ltr_all2 <- lm(I(100*log(meanco2))~year,data=g3)
  } else if (ratio=="logodds")
  {
    ltr <- lm(I(100*dv)~country+country:year,data=m[(m$country %in% "France" & m$year==1993)==F
                                                    & (m$country %in% "Denmark" & m$year==1994)==F,])
    ltr_all <- lm(I(100*dv)~country+year,data=m[(m$country %in% "France" & m$year==1993)==F
                                                & (m$country %in% "Denmark" & m$year==1994)==F,])
    ltr_all2 <- lm(I(100*meanco2)~year,data=g3)
  }
  print(screenreg(list(ltr,ltr_all,ltr_all2),digits=4))
  
  return(list(m=m,g=g,g2=g2,g3=g3,h2=h2,mycol=mycol,des=des
              ,ltr=ltr,ltr_all=ltr_all,ltr_all=ltr_all2
  ))
  
}

#---------------------------------#
# FIGURES ####
#---------------------------------#

# Country definition and colors

es="Spain"
sw="Sweden"
ca="Canada"
fr="France"
no="Norway"
dk="Denmark"
jp="Japan"
de="Germany"
cz="Czechia"
kr="South Korea"
hu="Hungary"
nl="Netherlands"

Canada <- c(brewer.pal(8, "Blues")[8],"1","21")
USA <- c(brewer.pal(8, "Blues")[6],"3","1")

Sweden <- c(brewer.pal(8, "BuGn")[8],"1","2")
Norway <- c(brewer.pal(8, "Dark2")[5],"1","24")
Denmark <- c(brewer.pal(8, "BuGn")[6],"1","25")

France <- c(brewer.pal(8, "RdPu")[8],"1","22")
Netherlands<-c(brewer.pal(8, "RdPu")[6],"1","7")
Germany <- c(brewer.pal(8, "RdPu")[4],"2","0")

Spain <- c(brewer.pal(8, "Set1")[1],"2","8")

Czechia <- c(brewer.pal(8, "BrBG")[1],"5","23")
Hungary <- c(brewer.pal(8, "BrBG")[2],"5","9")

Japan <- c(brewer.pal(8, "Oranges")[6],"2","4")
Korea <- c(brewer.pal(8, "Oranges")[5],"2","3")

ctrycolor <- data.frame(Canada,Czechia,Denmark,France,Germany,Hungary,Japan,Korea,Netherlands,Norway,Spain,Sweden)
colnames(ctrycolor)[which(names(ctrycolor) %in% "SouthKorea")] <- "South Korea"


## Figure 1. Top earners workplace isolation ####

#Data
all_e <- readRDS("Data/rds/exp_all_e.rds")
all_e5 <- all_e[all_e$fwage %in% 5,]
all_e45 <- all_e[all_e$fwage %in% 45,]
all_e1 <- all_e[all_e$fwage %in% 1,]

myctry <- names(table(all_e5$country))
mycol <- as.vector(na.omit(t(ctrycolor[1,which(names(ctrycolor) %in% myctry)])))
mylty <- as.vector(as.numeric(na.omit(t(ctrycolor[2,which(names(ctrycolor) %in% myctry)]))))
mypch <- as.vector(as.numeric(na.omit(t(ctrycolor[3,which(names(ctrycolor) %in% myctry)]))))

myorder <- c(1,NA,3,10,12,4,9,5,11,2,6,8,7)
myctry <- myctry[myorder]
mycol <- mycol[myorder]
mylty <- mylty[myorder]
mypch <- mypch[myorder]

#Figure
pdf(file="Figures/Figure_1_top_earner_isol.pdf",width=7,height=10,family="serif")
par(mfrow=c(2,1))
g_e55 <- mygraph_f(all_e5$mean_xF9910,"",
                   "Exposure (log-odds scale)","",ratio="OR",legpos="top",myyaxt="n",
                   my_mar=c(0.1,4,3,1),my_inset=c(0,-0.1),myxaxt="n",my_legend = "n")
axis(2,logodds100(c(11,13,15,17,19,21,22)),label = paste(c(11,13,15,17,19,21,22),"%"), las=2,cex.axis=0.8)
text(1991, logodds100(20), "A. Top 1%",cex = .8)

if (g_e55$ltr_all$coefficients["year"]>0){plus <- "+"
}  else {plus <- ""}
my_delta <- paste("\n",plus,round(g_e55$ltr_all$coefficients["year"],1),"% /year",sep="")

text(g_e55$des$mystart[g_e55$des$country %in% "Adj. Mean"],
     logodds100(-0.2+g_e55$des$first[g_e55$des$country %in% "Adj. Mean"]), 
     bquote (Delta : .(my_delta)),cex = .6,pos=1)
legend("top", 
       inset=c(0,-0.15),
       legend=c("Adj. Mean",myctry),
       pch =c(NA,mypch), 
       col=c("black",mycol),
       pt.bg=c("black",mycol),
       lty=c(1,mylty),
       lwd = c(3,rep(1,length(myctry))),
       text.col = c("black",mycol),
       cex=0.8,ncol=6,bg=NA,box.col="gray",box.lty=0)
g_e4545 <- mygraph_f(all_e45$mean_xF9010,"",
                     "Exposure (log-odds scale)","",ratio="OR",legpos = "top",myyaxt="n",
                     my_mar=c(2.5,4,0.1,1),my_inset=c(0,-0.05),my_legend = "n")
#my_ylim = c(logodds100(10.5),logodds100(40))
axis(2,logodds100(c(22:24,26:29,31:34,36:39,41:42)),label = paste(c(22:24,26:29,31:34,36:39,41:42),"%"), las=2,cex.axis=0.8)
text(1991, logodds100(39), "B. Top 10%",cex = .8)
if (g_e4545$ltr_all$coefficients["year"]>0){plus <- "+"
}  else {plus <- ""}
my_delta <- paste("\n",plus,round(g_e4545$ltr_all$coefficients["year"],1),"% /year",sep="")
text(g_e4545$des$mystart[g_e4545$des$country %in% "Adj. Mean"],
     logodds100(-0.2+g_e4545$des$first[g_e4545$des$country %in% "Adj. Mean"]), 
     bquote (Delta : .(my_delta)),cex = .6,pos=1)
dev.off()


## Figure 2. Top earners workplace exposure to bottom 25% ####
pdf(file="Figures/Figure_2_top_earn_exp_bott.pdf",width=7,height=10,family="serif")
par(mfrow=c(2,1))
g_e51 <- mygraph_f(all_e5$mean_xF0025,"",
                   "Exposure to bottom 25% (log-odds scale)","",ratio="OR",legpos = "top",myyaxt="n",
                   my_mar=c(0.1,4,3,1),my_inset=c(0,-0.1),myxaxt="n",my_ylim=c(logodds100(1.8),logodds100(16)),
                   my_legend = "n")
axis(2,logodds100(c(11,13,15)),label = paste(c(11,13,15),"%"), las=2,cex.axis=0.8)
text(1991, logodds100(2), "A. Top 1%",cex = .8)

if (g_e51$ltr_all$coefficients["year"]>0){plus <- "+"
}  else {plus <- ""}
my_delta <- paste("\n",plus,round(g_e51$ltr_all$coefficients["year"],1),"% /year",sep="")
text(g_e51$des$mystart[g_e51$des$country %in% "Adj. Mean"],
     logodds100(-0.2+g_e51$des$first[g_e51$des$country %in% "Adj. Mean"]), 
     bquote (Delta : .(my_delta)),cex = .6,pos=1)
legend("top", 
       inset=c(0,-0.15),
       legend=c("Adj. Mean",myctry),pch =c(NA,mypch), 
       col=c("black",mycol),
       pt.bg=c("black",mycol),cex=0.8,ncol=6,lty=c(1,mylty),lwd = c(3,rep(1,length(myctry))),
       text.col = c("black",mycol),bg=NA,box.col="gray",box.lty=0)
g_e451 <- mygraph_f(all_e45$mean_xF0025,"",
                    "Exposure to bottom 25%  (log-odds scale)","",ratio="OR",legpos = "topright",myyaxt="n",
                    my_mar=c(2.5,4,0.1,1),my_inset=c(0,-0.05),my_legend = "n",my_ylim=c(logodds100(1.8),logodds100(16)))
axis(2,logodds100(c(11,13,15)),label = paste(c(11,13,15),"%"), las=2,cex.axis=0.8
)

if (g_e451$ltr_all$coefficients["year"]>0){plus <- "+"
}  else {plus <- ""}
my_delta <- paste("\n",plus,round(g_e451$ltr_all$coefficients["year"],1),"% /year",sep="")
text(1991, logodds100(2), "B. Top 10%",cex = .8)
text(g_e451$des$mystart[g_e451$des$country %in% "Adj. Mean"],
     logodds100(-0.2+g_e451$des$first[g_e451$des$country %in% "Adj. Mean"]), 
     bquote (Delta : .(my_delta)),cex = .6,pos=1)
dev.off()


## Figure 3. Bottom 25%  workplace isolation ####
pdf(file="Figures/Figure_3_bottom_isol.pdf",width=7,height=6,family="serif")
g_e11 <- mygraph_f(all_e1$mean_xF0025,"","Exposure (% - log-odds scale)","",ratio="OR",legpos="top",myyaxt="n",
                   my_mar=c(2.5,4,3,1),my_inset=c(0,-0.1),my_legend = "n")                
#axis(2,logodds100(c(36:39,41:44,46:49,51:54,56:59,61:64)),label = paste(c(36:39,41:44,46:49,51:54,56:59,61:64),"%"), las=2,cex.axis=0.8)
#text(1990, logodds100(65), "Bottom 25%",cex = .8)

if (g_e11$ltr_all$coefficients["year"]>0){plus <- "+"
}  else {plus <- ""}
my_delta <- paste("\n",plus,round(g_e11$ltr_all$coefficients["year"],1),"% /year",sep="")
text(g_e11$des$mystart[g_e11$des$country %in% "Adj. Mean"],
     logodds100(-0.2+g_e11$des$first[g_e11$des$country %in% "Adj. Mean"]), 
     bquote (Delta : .(my_delta)),cex = .6,pos=1)
legend("top", 
       inset=c(0,-0.12),
       legend=c("Adj. Mean",myctry),pch =c(NA,mypch), 
       col=c("black",mycol),pt.bg=c("black",mycol),
       cex=0.8,ncol=6,lty=c(1,mylty),lwd = c(3,rep(1,length(myctry))),
       text.col = c("black",mycol),bg=NA,box.col="gray",box.lty=0)
dev.off()


## Figure 4. Decile exposure to one another  ####
#---------------------------##

# Data and preliminary treatments
d_all <- readRDS("Data/rds/exp_d_all.rds")


dd0<-felm(I(100*logodds(mean_xd0))~country:factor(d):year|paste(factor(d),country)|0|year,
          data=d_all[(d_all$country %in% "France" & d_all$year==1993)==F
                     & (d_all$country %in% "Denmark" & d_all$year==1994)==F,])
screenreg(dd0,digits=2)


dd1<-felm(I(100*logodds(mean_xd1))~country:factor(d):year|paste(factor(d),country)|0|year,
          data=d_all[(d_all$country %in% "France" & d_all$year==1993)==F
                     & (d_all$country %in% "Denmark" & d_all$year==1994)==F,])
screenreg(dd1,digits=2)


dd2 <- felm(I(100*logodds(mean_xd2))~country:factor(d):year|paste(factor(d),country)|0|year,
            data=d_all[(d_all$country %in% "France" & d_all$year==1993)==F
                       & (d_all$country %in% "Denmark" & d_all$year==1994)==F,])
screenreg(dd2,digits=2)

dd3 <- felm(I(100*logodds(mean_xd3))~country:factor(d):year|paste(factor(d),country)|0|year,
            data=d_all[(d_all$country %in% "France" & d_all$year==1993)==F
                       & (d_all$country %in% "Denmark" & d_all$year==1994)==F,])
screenreg(dd3,digits=2)

dd4 <- felm(I(100*logodds(mean_xd4))~country:factor(d):year|paste(factor(d),country)|0|year,
            data=d_all[(d_all$country %in% "France" & d_all$year==1993)==F
                       & (d_all$country %in% "Denmark" & d_all$year==1994)==F,])
screenreg(dd4,digits=2)

dd5 <- felm(I(100*logodds(mean_xd5))~country:factor(d):year|paste(factor(d),country)|0|year,
            data=d_all[(d_all$country %in% "France" & d_all$year==1993)==F
                       & (d_all$country %in% "Denmark" & d_all$year==1994)==F,])
screenreg(dd5,digits=2)

dd6 <- felm(I(100*logodds(mean_xd6))~country:factor(d):year|paste(factor(d),country)|0|year,
            data=d_all[(d_all$country %in% "France" & d_all$year==1993)==F
                       & (d_all$country %in% "Denmark" & d_all$year==1994)==F,])
screenreg(dd6,digits=2)

dd7 <- felm(I(100*logodds(mean_xd7))~country:factor(d):year|paste(factor(d),country)|0|year,
            data=d_all[(d_all$country %in% "France" & d_all$year==1993)==F
                       & (d_all$country %in% "Denmark" & d_all$year==1994)==F,])
screenreg(dd7,digits=2)

dd8 <- felm(I(100*logodds(mean_xd8))~country:factor(d):year|paste(factor(d),country)|0|year,
            data=d_all[(d_all$country %in% "France" & d_all$year==1993)==F
                       & (d_all$country %in% "Denmark" & d_all$year==1994)==F,])
screenreg(dd8,digits=2)

dd9 <- felm(I(100*logodds(mean_xd9))~country:factor(d):year|paste(factor(d),country)|0|year,
            data=d_all[(d_all$country %in% "France" & d_all$year==1993)==F
                       & (d_all$country %in% "Denmark" & d_all$year==1994)==F,])
screenreg(dd9,digits=2)
str(dd9)

dd_all0<-felm(I(100*logodds(mean_xd0))~factor(d):year|paste(factor(d),country)|0|year,
              data=d_all[(d_all$country %in% "France" & d_all$year==1993)==F
                         & (d_all$country %in% "Denmark" & d_all$year==1994)==F,])
screenreg(dd_all0,digits=2)


dd_all1<-felm(I(100*logodds(mean_xd1))~factor(d):year|paste(factor(d),country)|0|year,
              data=d_all[(d_all$country %in% "France" & d_all$year==1993)==F
                         & (d_all$country %in% "Denmark" & d_all$year==1994)==F,])
screenreg(dd_all1,digits=2)


dd_all2 <- felm(I(100*logodds(mean_xd2))~factor(d):year|paste(factor(d),country)|0|year,
                data=d_all[(d_all$country %in% "France" & d_all$year==1993)==F
                           & (d_all$country %in% "Denmark" & d_all$year==1994)==F,])
screenreg(dd_all2,digits=2)

dd_all3 <- felm(I(100*logodds(mean_xd3))~factor(d):year|paste(factor(d),country)|0|year,
                data=d_all[(d_all$country %in% "France" & d_all$year==1993)==F
                           & (d_all$country %in% "Denmark" & d_all$year==1994)==F,])
screenreg(dd_all3,digits=2)

dd_all4 <- felm(I(100*logodds(mean_xd4))~factor(d):year|paste(factor(d),country)|0|year,
                data=d_all[(d_all$country %in% "France" & d_all$year==1993)==F
                           & (d_all$country %in% "Denmark" & d_all$year==1994)==F,])
screenreg(dd_all4,digits=2)

dd_all5 <- felm(I(100*logodds(mean_xd5))~factor(d):year|paste(factor(d),country)|0|year,
                data=d_all[(d_all$country %in% "France" & d_all$year==1993)==F
                           & (d_all$country %in% "Denmark" & d_all$year==1994)==F,])
screenreg(dd_all5,digits=2)

dd_all6 <- felm(I(100*logodds(mean_xd6))~factor(d):year|paste(factor(d),country)|0|year,
                data=d_all[(d_all$country %in% "France" & d_all$year==1993)==F
                           & (d_all$country %in% "Denmark" & d_all$year==1994)==F,])
screenreg(dd_all6,digits=2)

dd_all7 <- felm(I(100*logodds(mean_xd7))~factor(d):year|paste(factor(d),country)|0|year,
                data=d_all[(d_all$country %in% "France" & d_all$year==1993)==F
                           & (d_all$country %in% "Denmark" & d_all$year==1994)==F,])
screenreg(dd_all7,digits=2)

dd_all8 <- felm(I(100*logodds(mean_xd8))~factor(d):year|paste(factor(d),country)|0|year,
                data=d_all[(d_all$country %in% "France" & d_all$year==1993)==F
                           & (d_all$country %in% "Denmark" & d_all$year==1994)==F,])
screenreg(dd_all8,digits=2)

dd_all9 <- felm(I(100*logodds(mean_xd9))~factor(d):year|paste(factor(d),country)|0|year,
                data=d_all[(d_all$country %in% "France" & d_all$year==1993)==F
                           & (d_all$country %in% "Denmark" & d_all$year==1994)==F,])
screenreg(dd_all9,digits=2)
str(dd_all9)


dd <- data.frame(dd0$coefficients,
                 dd1$coefficients,
                 dd2$coefficients,
                 dd3$coefficients,
                 dd4$coefficients,
                 dd5$coefficients,
                 dd6$coefficients,
                 dd7$coefficients,
                 dd8$coefficients,
                 dd9$coefficients)
colnames(dd) <- c("xd0","xd1","xd2","xd3","xd4","xd5","xd6","xd7","xd8","xd9")
dd$country <- gsub("country","",substr(rownames(dd),1,regexpr(":",rownames(dd))-1))
dd$d <- gsub(":year","",substr(rownames(dd),regexpr(":",rownames(dd))+1+9,nchar(rownames(dd))))
rownames(dd) <- NULL

dd_all <- data.frame(dd_all0$coefficients,
                     dd_all1$coefficients,
                     dd_all2$coefficients,
                     dd_all3$coefficients,
                     dd_all4$coefficients,
                     dd_all5$coefficients,
                     dd_all6$coefficients,
                     dd_all7$coefficients,
                     dd_all8$coefficients,
                     dd_all9$coefficients)
colnames(dd_all) <- c("xd0","xd1","xd2","xd3","xd4","xd5","xd6","xd7","xd8","xd9")
dd_all$d <- substr(rownames(dd_all),10,10)
table(dd_all$d)
dd_all$country <- "All countries"
rownames(dd_all) <- NULL
dd <- rbind(dd,dd_all)

#Figure
dec_evol <- function(ctry,start,end,scale=2) {    
  
  mygr_rate <- dd[dd$country==ctry,]
  
  pdf_name <- paste("Figures/Figure_4_",ctry,".pdf",sep="")
  pdf(file=pdf_name,paper="a4r",width=10,height=7,family="serif")
  par(xpd=FALSE, mar=c(2.5,2.5,1,1))
  matplot(10*c(0:9)+10,cbind(diag(as.matrix(mygr_rate[,c(1:10)])),
                             mygr_rate[,c(1:10)]),
          pch=c(1,0,2,5,6,3,4,25,23,24,22),
          cex=c(3,rep(1,10)),
          col=c("black",spectral[-6]),
          bg=c("black",spectral[-6]),
          xlab=ctry,
          las = 2, 
          axes = FALSE,
          cex.lab=0.8,
          cex.axis=0.8,
          xlim=c(0,100),#ylim=c(0.001,0.5),
          ylab=paste(start,"-",end," yearly trend in rate of evolution of deciles' exposure to one another (%)", sep="")
  )
  axis(side = 1, at = c(0:9)*10+10 ,label =c("D1","D2","D3","D4","D5","D6","D7","D8","D9","D10")
  )
  axis(side = 1, at = c(-5,10),label=c("",""))
  axis(side = 2, at = c(-20:20)/scale)
  title(ctry)
  matlines(10*c(0:9)+10,mygr_rate[,c(1:10)],lwd = 1,col=spectral[-6],lty=1)
  matlines(-5:110,c(rep(0,116)),lwd = 1,lty=5,col="gray")
  legend("bottomleft",bg="transparent",legend=c("D1","D2","D3","D4","D5","D6","D7","D8","D9","D10"),
         text.col = spectral[-6],box.lwd=-1, lty=1,
         pch=c(0,2,5,6,3,4,25,23,24,22),
         pt.bg=c(spectral[-6]),
         ncol=1,cex=1,col=spectral[-6] #,inset=c(0,-0.1)
  )
  dev.off()
  

}

dec_evol(ctry="All countries",start=1990,end=2019) 

## Figure 5 French events####

fr_ev <- readRDS("Data/rds/exp_fr_ev.rds")
dd_out <- fr_ev[fr_ev$event=="out",]

# Plot 
plot_out <- ggplot(dd_out, aes(x=time, y=b,color=ctrl)) +        # ggplot2 plot with confidence intervals
  theme_bw()+
  theme(panel.grid.major = element_line(color = "white"))+
  geom_hline(yintercept=0,color="black")+
  ylab("Effect of 10 percent workforce outsourcing on top 10% isolation (in percentage point)") +
  geom_vline(xintercept=-0.5,color="black",linetype="dashed",size=1)+
  # geom_vline(xintercept=2011.5,color="gray5",linetype="dashed")+
  geom_point(aes(shape=ctrl),size=3) +
  geom_line() +
  
  geom_errorbar(aes(ymin = l, ymax = u),width=0.5)+
  scale_x_continuous(breaks = round(seq(min(floor(dd_out$time)), max(dd_out$time), by = 1),1))+
  # annotate("text", x=2010, y=-2.0, label= "Offshoring event
  #         between 2009 & 2011",cex=3) +
  # annotate("text", x=2006.5, y=5, label= "Pre-event difference",col="red") +
  # annotate("text", x=2013.5, y=-2, label= "Post-event difference",col="blue") +
  scale_color_manual(values=c("darkblue", "darkorange"))+
  theme(legend.position="bottom")+ 
  theme(legend.title = element_blank())+
  ylab("Top 10% isolation (in % point)") +
  ggtitle("Outsourcing (10% of workforce)")+
  theme(legend.position="none") +
  # theme(axis.title.y = element_blank()) +
  theme(axis.title.x = element_blank())

# dev.off()
plot_out

dd_lay <- fr_ev[fr_ev$event=="lay",]
# Plot 
plot_lay <- ggplot(dd_lay, aes(x=time, y=b,color=ctrl)) +        # ggplot2 plot with confidence intervals
  theme_bw()+
  theme(panel.grid.major = element_line(color = "white"))+
  geom_hline(yintercept=0,color="black")+
  geom_vline(xintercept=-0.5,color="black",linetype="dashed",size=1)+
  # geom_vline(xintercept=2011.5,color="gray5",linetype="dashed")+
  geom_point(aes(shape=ctrl),size=3) +
  geom_line() +
  geom_errorbar(aes(ymin = l, ymax = u),width=0.5)+
  scale_x_continuous(breaks = round(seq(min(floor(dd_lay$time)), max(dd_lay$time), by = 1),1))+
  # annotate("text", x=2010, y=-2.0, label= "Offshoring event
  #         between 2009 & 2011",cex=3) +
  # annotate("text", x=2006.5, y=5, label= "Pre-event difference",col="red") +
  # annotate("text", x=2013.5, y=-2, label= "Post-event difference",col="blue") +
  theme(legend.position="bottom")+
  theme(legend.title = element_blank())+
  theme(legend.position="none") +
  ylab("Top 10% isolation (in % point)") +
  # theme(axis.title.y = element_blank()) +
  # theme(axis.title.x = element_blank()) +
  scale_color_manual(values=c("darkblue", "darkorange"))+
  ggtitle("Layoffs (10% of workforce)") 
dev.off()
plot_lay


# Offshoring
dd_off <- fr_ev[fr_ev$event=="off",]
plot_off <- ggplot(dd_off, aes(x=time, y=b,color=ctrl)) +        # ggplot2 plot with confidence intervals
  theme_bw()+
  theme(panel.grid.major = element_line(color = "white"))+
  geom_hline(yintercept=0,color="black")+
  # ylab("Effect of 10 percent workforce offshoring on top 10% isolation (in percentage point)") +
  geom_vline(xintercept=2008.5,color="black",linetype="dashed",size=1)+
  geom_vline(xintercept=2011.5,color="gray5",linetype="dashed")+
  geom_point(aes(shape=ctrl),size=3) +
  geom_line() +
  geom_errorbar(aes(ymin = l, ymax = u),width=0.5)+
  scale_x_continuous(breaks = round(seq(min(floor(dd_off$time)), max(dd_off$time), by = 1),1))+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  annotate("text", x=2010, y=-2.0, label= "Event between
  2009 & 2011",cex=3) +
  # annotate("text", x=2006.5, y=5, label= "Pre-event difference",col="red") +
  # annotate("text", x=2013.5, y=-2, label= "Post-event difference",col="blue") +
  theme(legend.position="bottom")+ 
  theme(legend.title = element_blank())+
  scale_color_manual(values=c("darkblue", "darkorange"))+
  ylab("Top 10% isolation (in % point)") +
  ggtitle("Offshoring (10% of workforce)")+
  theme(legend.position="none") +
  theme(axis.title.y = element_blank()) +
  theme(axis.title.x = element_blank())

dev.off()
plot_off


# Subcontracting
dd_sub <- fr_ev[fr_ev$event=="sub",]
plot_sub <- ggplot(dd_sub, aes(x=time, y=b,color=ctrl)) +        # ggplot2 plot with confidence intervals
  theme_bw()+
  theme(panel.grid.major = element_line(color = "white"))+
  geom_hline(yintercept=0,color="black")+
  # ylab("Effect of 10 percent workforce offshoring on top 10% isolation (in percentage point)") +
  geom_vline(xintercept=2011.5,color="black",linetype="dashed",size=1)+
  geom_vline(xintercept=2016.5,color="gray5",linetype="dashed")+
  geom_point(aes(shape=ctrl),size=3) +
  geom_line() +
  geom_errorbar(aes(ymin = l, ymax = u),width=0.5)+
  scale_x_continuous(breaks = round(seq(min(floor(dd_sub$time)), max(dd_sub$time), by = 1),1))+
  annotate("text", x=2014, y=-1, label= "2011-2017 change
  in subcontracting",cex=3) +
  # annotate("text", x=2006.5, y=5, label= "Pre-event difference",col="red") +
  # annotate("text", x=2013.5, y=-2, label= "Post-event difference",col="blue") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  theme(legend.position="bottom")+ 
  theme(legend.title = element_blank())+
  ylab("Top 10% isolation (in % point)") +
  ggtitle("Subcontracting (10% of turnover)")+
  scale_color_manual(values=c("darkblue", "darkorange"))+
  # theme(legend.position="none") +
  theme(axis.title.y = element_blank()) +
  theme(axis.title.x = element_blank())
dev.off()
plot_sub


grid.arrange(plot_out, plot_off, plot_lay, plot_sub, ncol=2)

pdf("Figures/Figure5_France_wgt_event.pdf")
grid.arrange(plot_out, plot_off, plot_lay, plot_sub, ncol=2)
dev.off()



# Supplementary figures ####
## Figure S3.1. Top earners isolation with confidence intervals ####


  # Prepatory work 
  e55_l <- mygraph2(all_e5$ic9910l95,"Top 1% establishment isolation",
                    "Exposure (% - log-odds scale)","",ratio="OR",legpos="bottomleft",myyaxt="n")
  e55_u <- mygraph2(all_e5$ic9910u95,"Top 1% establishment isolation",
                    "Exposure (% - log-odds scale)","",ratio="OR",legpos="bottomleft",myyaxt="n")
  
  d_all$icd9l95 <- d_all$mean_xd9-1.96*d_all$sd_xd9/(d_all$N_xd9)**0.5
  d_all$icd9u95 <- d_all$mean_xd9+1.96*d_all$sd_xd9/(d_all$N_xd9)**0.5
  
  d9 <- d_all[d_all$d %in% 9,]
  
  d9_l <- mygraph2(d9$icd9l95,"Top 10% establishment isolation",
                   "Exposure (% - log-odds scale)","",ratio="OR",legpos="bottomleft",myyaxt="n")
  d9_u <- mygraph2(d9$icd9u95,"Top 10% establishment isolation",
                   "Exposure (% - log-odds scale)","",ratio="OR",legpos="bottomleft",myyaxt="n")

#figure
pdf(file="Figures/Figure_S3_1.pdf",width=7,height=10,family="serif")
par(mfrow=c(2,1))
z <- mygraph_f(all_e5$mean_xF9910,"",
               "Exposure (log-odds scale)","",ratio="OR",legpos="top",myyaxt="n",myxaxt="n",
               my_mar=c(0.1,4,3,1),
               my_inset=c(0,-0.1),
               my_legend = "n")
matlines(e55_l$g2$year,logodds100(cbind(e55_l$g2[,-1],e55_u$g2[,-1])),col=rep(z$mycol,2),lty=2)
axis(2,logodds100(c(11,13,15)),label = paste(c(11,13,15),"%"), las=2,cex.axis=0.8)
text(1991, logodds100(6), "A. Top 1%",cex = .8)
legend("top", 
       inset=c(0,-0.15),
       legend=c("Adj. Mean",myctry),pch =c(NA,mypch), 
       col=c("black",mycol),pt.bg=c("black",mycol),
       cex=0.8,ncol=6,lty=c(1,mylty),lwd = c(3,rep(1,length(myctry))),
       text.col = c("black",mycol),bg=NA,box.col="gray",box.lty=0)
g_d9 <- mygraph_f(d9$mean_xd9,"","Exposure (log-odds scale)","",
                  ratio="OR",legpos="topright",myyaxt="n",
                  my_legend = "n",
                  my_mar=c(2.5,4,0.1,1),my_inset=c(0,-0.2))
matlines(d9_l$g2$year,logodds100(cbind(d9_l$g2[,-1],d9_u$g2[,-1])),col=rep(z$mycol,2),lty=2)
text(1991, logodds100(39), "B. Top 10%",cex = .8)
dev.off()

## Figure S3.2. Top earners workplace exposure to mid quartiles ####

pdf(file="Figures/Figure_S3_2.pdf",width=7,height=10,family="serif")
par(mfrow=c(2,1))
f1 <- mygraph_f(all_e5$mean_xF2575,"",
                "Exposure (log-odds scale)","",ratio="OR",legpos = "top",myyaxt="n",
                my_mar=c(0.1,4,3,1),my_inset=c(0,-0.1),myxaxt="n",my_legend = "n"
                #my_ylim=c(logodds100(1.8),logodds100(16))
)
axis(2,logodds100(c(11,13,15)),label = paste(c(11,13,15),"%"), las=2,cex.axis=0.8)
text(1991, logodds100(20), "A. Top 1%",cex = .8)
if (f1$ltr_all$coefficients["year"]>0){plus <- "+"
}  else {plus <- ""}
my_delta <- paste("\n",plus,round(f1$ltr_all$coefficients["year"],1),"% /year",sep="")
text(f1$des$mystart[f1$des$country %in% "Adj. Mean"],
     logodds100(-0.2+f1$des$first[f1$des$country %in% "Adj. Mean"]), 
     bquote (Delta : .(my_delta)),cex = .6,pos=1)
legend("top", 
       inset=c(0,-0.15),
       legend=c("Adj. Mean",myctry),pch =c(NA,mypch), 
       col=c("black",mycol),pt.bg=c("black",mycol),cex=0.8,ncol=6,lty=c(1,mylty),lwd = c(3,rep(1,length(myctry))),
       text.col = c("black",mycol),bg=NA,box.col="gray",box.lty=0)
f2 <- mygraph_f(all_e45$mean_xF2575,"",
                "Exposure (log-odds scale)","",ratio="OR",legpos = "topright",myyaxt="n",
                my_mar=c(2.5,4,0.1,1),my_inset=c(0,-0.05),my_legend = "n"
                #my_ylim=c(logodds100(1.8),logodds100(16))
)
axis(2,logodds100(c(11,13,15)),label = paste(c(11,13,15),"%"), las=2,cex.axis=0.8
)
text(1991, logodds100(30), "B. Top 10%",cex = .8)
if (f2$ltr_all$coefficients["year"]>0){plus <- "+"
}  else {plus <- ""}
my_delta <- paste("\n",plus,round(f2$ltr_all$coefficients["year"],1),"% /year",sep="")
text(f2$des$mystart[f2$des$country %in% "Adj. Mean"],
     logodds100(-0.2+f2$des$first[f2$des$country %in% "Adj. Mean"]), 
     bquote (Delta : .(my_delta)),cex = .6,pos=1)
dev.off()

## Figure S3.3 Country decile exposures ####
dec_evol_s <- function(ctry,start,end,scale=2) {    
  
  mygr_rate <- dd[dd$country==ctry,]
  
  pdf_name <- paste("Figures/Figure_S3_3_",ctry,".pdf",sep="")
  pdf(file=pdf_name,paper="a4r",width=10,height=7,family="serif")
  par(xpd=FALSE, mar=c(2.5,2.5,1,1))
  matplot(10*c(0:9)+10,cbind(diag(as.matrix(mygr_rate[,c(1:10)])),
                             mygr_rate[,c(1:10)]),
          pch=c(1,0,2,5,6,3,4,25,23,24,22),
          cex=c(3,rep(1,10)),
          col=c("black",spectral[-6]),
          bg=c("black",spectral[-6]),
          xlab=ctry,
          las = 2, 
          axes = FALSE,
          cex.lab=0.8,
          cex.axis=0.8,
          xlim=c(0,100),#ylim=c(0.001,0.5),
          ylab=paste(start,"-",end," yearly trend in rate of evolution of deciles' exposure to one another (%)", sep="")
  )
  axis(side = 1, at = c(0:9)*10+10 ,label =c("D1","D2","D3","D4","D5","D6","D7","D8","D9","D10")
  )
  axis(side = 1, at = c(-5,10),label=c("",""))
  axis(side = 2, at = c(-20:20)/scale)
  title(ctry)
  matlines(10*c(0:9)+10,mygr_rate[,c(1:10)],lwd = 1,col=spectral[-6],lty=1)
  matlines(-5:110,c(rep(0,116)),lwd = 1,lty=5,col="gray")
  legend("bottomleft",bg="transparent",legend=c("D1","D2","D3","D4","D5","D6","D7","D8","D9","D10"),
         text.col = spectral[-6],box.lwd=-1, lty=1,
         pch=c(0,2,5,6,3,4,25,23,24,22),
         pt.bg=c(spectral[-6]),
         ncol=1,cex=1,col=spectral[-6] #,inset=c(0,-0.1)
  )
  dev.off()
  
  pdf_name <- paste("Figures//Figure_S3_3_",ctry,".pdf",sep="")
  pdf(file=pdf_name,width=10,height=7,family="serif")
  
  par(xpd=FALSE, mar=c(2.5,2.5,1,1))
  matplot(10*c(0:9)+10,cbind(diag(as.matrix(mygr_rate[,c(1:10)])),
                             mygr_rate[,c(1:10)]),
          pch=c(1,0,2,5,6,3,4,25,23,24,22),
          cex=c(3,rep(1,10)),
          col=c("black",spectral[-6]),
          bg=c("black",spectral[-6]),
          xlab=ctry,
          las = 2, 
          axes = FALSE,
          cex.lab=0.8,
          cex.axis=0.8,
          xlim=c(0,100),#ylim=c(0.001,0.5),
          ylab=paste(start,"-",end," yearly trend in rate of evolution of deciles' exposure to one another (%)", sep="")
  )
  axis(side = 1, at = c(0:9)*10+10 ,label =c("D1","D2","D3","D4","D5","D6","D7","D8","D9","D10")
  )
  axis(side = 1, at = c(-5,10),label=c("",""))
  axis(side = 2, at = c(-20:20)/scale)
  title(ctry)
  matlines(10*c(0:9)+10,mygr_rate[,c(1:10)],lwd = 1,col=spectral[-6],lty=1)
  matlines(-5:110,c(rep(0,116)),lwd = 1,lty=5,col="gray")
  legend("bottomleft",bg="transparent",legend=c("D1","D2","D3","D4","D5","D6","D7","D8","D9","D10"),
         text.col = spectral[-6],box.lwd=-1, lty=1,
         pch=c(0,2,5,6,3,4,25,23,24,22),
         pt.bg=c(spectral[-6]),
         ncol=1,cex=1,col=spectral[-6] #,inset=c(0,-0.1)
  )
  dev.off()
  
  

}

dec_evol_s(ctry="Canada",start=1990,end=2019,scale=4)

dec_evol_s(ctry="Czechia",start=2002,end=2016)

dec_evol_s(ctry="Denmark",start=1995,end=2018)

dec_evol_s(ctry="France",start=1994,end=2019)

dec_evol_s(ctry="Germany",start=1999,end=2015)

dec_evol_s(ctry="Hungary",start=2003,end=2017,scale=1)

dec_evol_s(ctry="Japan",start=1990,end=2013)

dec_evol_s(ctry="Korea",start=1990,end=2012)

dec_evol_s(ctry="Norway",start=1995,end=2018)

dec_evol_s(ctry="Netherlands",start=2006,end=2018)

dec_evol_s(ctry="Spain",start=2006,end=2018)



## Figure S3.4. Top earners share ####
top1share <- 0.01*all_e[all_e$fwage %in% 5,"mean_wage"]/all_e[all_e$fwage %in% NA,"mean_wage"]
top10share <- 0.1*(0.9*all_e[all_e$fwage %in% 4,"mean_wage"]+0.1*all_e[all_e$fwage %in% 5,"mean_wage"])/all_e[all_e$fwage %in% NA,"mean_wage"]

all_eshare <- data.frame(top1share,top10share,all_e[all_e$fwage %in% NA,c("mean_wage","country","year")])


pdf(file="Figures/Figure_S3_4.pdf",width=7,height=10,family="serif")
par(mfrow=c(2,1))
f1 <- mygraph_f(all_eshare$top1share,"",
                "Earnings share (log-odds scale)","",ratio="OR",legpos = "top",myyaxt="n",
                my_mar=c(0.1,4,3,1),my_inset=c(0,-0.1),myxaxt="n",my_legend = "n"
                #my_ylim=c(logodds100(1.8),logodds100(16))
)
axis(2,logodds100(c(11)),label = paste(c(11),"%"), las=2,cex.axis=0.8)
text(1991, logodds100(9), "A. Top 1%",cex = .8)
if (f1$ltr_all$coefficients["year"]>0){plus <- "+"
}  else {plus <- ""}
my_delta <- paste("\n",plus,round(f1$ltr_all$coefficients["year"],1),"% /year",sep="")
text(f1$des$mystart[f1$des$country %in% "Adj. Mean"],
     logodds100(-0.2+f1$des$first[f1$des$country %in% "Adj. Mean"]), 
     bquote (Delta : .(my_delta)),cex = .6,pos=1)
legend("top", 
       inset=c(0,-0.15),
       legend=c("Adj. Mean",myctry),pch =c(NA,mypch), 
       col=c("black",mycol),
       pt.bg=c("black",mycol),cex=0.8,ncol=6,lty=c(1,mylty),lwd = c(3,rep(1,length(myctry))),
       text.col = c("black",mycol),bg=NA,box.col="gray",box.lty=0)
f2 <- mygraph_f(all_eshare$top10share,"",
                "Earnings share (log-odds scale)","",ratio="OR",legpos = "topright",myyaxt="n",
                my_mar=c(2.5,4,0.1,1),my_inset=c(0,-0.05),my_legend = "n"
                #my_ylim=c(logodds100(1.8),logodds100(16))
)
axis(2,logodds100(c(21:24,26:29,31:34)),label = paste(c(21:24,26:29,31:34),"%"), las=2,cex.axis=0.8
)
text(1991, logodds100(30), "B. Top 10%",cex = .8)
if (f2$ltr_all$coefficients["year"]>0){plus <- "+"
}  else {plus <- ""}
my_delta <- paste("\n",plus,round(f2$ltr_all$coefficients["year"],1),"% /year",sep="")
text(f2$des$mystart[f2$des$country %in% "Adj. Mean"],
     logodds100(-0.2+f2$des$first[f2$des$country %in% "Adj. Mean"]), 
     bquote (Delta : .(my_delta)),cex = .6,pos=1)
dev.off()

## Figure S3.5. Top 1% & 10%  residential isolation ####
all_l <- readRDS("Data/rds/exp_all_l.rds")
all_l5 <- all_l[all_l$fwage==5,]
all_l45 <- all_l[all_l$fwage==45,]
  
pdf(file="Figures/Figure_S3_5.pdf",width=7,height=10,family="serif")
par(mfrow=c(2,1))
g_l55 <- mygraph_f(all_l5$mean_xF9910,"","Exposure (% - log-odds scale)","",
                   ratio="OR",legpos="top",myyaxt="n",myxaxt="n",
                   my_mar=c(0.1,4,3,1),my_inset=c(0,-0.1),my_legend = "n")                
#axis(2,logodds100(c(36:39,41:44,46:49,51:54,56:59,61:64)),label = paste(c(36:39,41:44,46:49,51:54,56:59,61:64),"%"), las=2,cex.axis=0.8)
#text(1990, logodds100(65), "Bottom 25%",cex = .8)
text(1991, logodds100(4), "A. Top 1%",cex = .8)
legend("top", 
       inset=c(0,-0.13),
       legend=c("Adj. Mean",myctry[-12]),pch =c(NA,mypch[-12]), 
       col=c("black",mycol[-12]),pt.bg=c("black",mycol[-12]),cex=0.8,ncol=6,lty=c(1,mylty[-12]),lwd = c(3,rep(1,length(myctry)[-12])),
       text.col = c("black",mycol[-12]),bg=NA,box.col="gray",box.lty=0)
text(g_l55$des$mystart[g_l55$des$country %in% "Adj. Mean"],
     logodds100(-0.05+g_l55$des$first[g_l55$des$country %in% "Adj. Mean"]), 
     bquote (Delta : .(g_l55$mydelta)),cex = .6,pos=1)

g_l4545 <- mygraph_f(all_l45$mean_xF9010,"","Exposure (% - log-odds scale)","",
                     ratio="OR",myyaxt = "n",
                     my_mar=c(2.5,4,0.1,1),my_inset=c(0,-0.1),my_legend = "n")
text(g_l4545$des$mystart[g_l4545$des$country %in% "Adj. Mean"],
     logodds100(-0.05+g_l4545$des$first[g_l4545$des$country %in% "Adj. Mean"]), 
     bquote (Delta : .(g_l4545$mydelta)),cex = .6,pos=1)
text(1991, logodds100(17), "B. Top 10%",cex = .8)
dev.off()


## Figure S5.1 Migrants and natives absolute exposure to migrants ####

# Data import
seg_mig <- readRDS("Data/rds/exp_seg_mig.rds")

seg_mig1 <- seg_mig[seg_mig$migrant %in% 1,]
seg_mig0 <- seg_mig[seg_mig$migrant %in% 0,]

seg_mig1$omigrant <- oddsratio(seg_mig1$mean_xmigrant,seg_mig0$mean_xmigrant)
seg_mig1$lomigrant <- log(seg_mig1$omigrant)

#Figure
pdf(file="Figures/Figure_S5_1.pdf",width=7,height=10,family="serif")
par(mfrow=c(2,1))
f1 <- mygraph_f(seg_mig1$mean_xmigrant,"",
                "Exposure (logodds scale)","",
                ratio="OR",legpos="top",
                myyaxt="n",
                my_mar=c(0.1,4,3,1),my_inset=c(0,-0.1),myxaxt="n",my_legend = "n")
text(1992, logodds100(10), "A. Migrants",cex = .8)
if (f1$ltr_all$coefficients["year"]>0){plus <- "+"
}  else {plus <- ""}
my_delta <- paste("\n",plus,round(f1$ltr_all$coefficients["year"],1),"% /year",sep="")
text(f1$des$mystart[f1$des$country %in% "Adj. Mean"],
     logodds100(-0.2+f1$des$first[f1$des$country %in% "Adj. Mean"]), 
     bquote (Delta : .(my_delta)),cex = .6,pos=1)
legend("top", 
       inset=c(0,-0.15),
       legend=c("Adj. Mean",myctry[-c(10:13)]),
       pch =c(NA,mypch[-c(10:13)]), 
       col=c("black",mycol[-c(10:13)]),
       pt.bg=c("black",mycol[-c(10:13)]),
       lty=c(1,mylty[-c(10:13)]),
       lwd = c(3,rep(1,length(myctry[-c(10:13)]))),
       text.col = c("black",mycol[-c(10:13)]),
       cex=0.8,ncol=4,
       bg=NA,box.col="gray",box.lty=0)
f2 <- mygraph_f(seg_mig0$mean_xmigrant,"",
                "Exposure (logodds scale)","",
                ratio="OR",legpos="top",
                myyaxt="n",my_legend = "n",
                my_mar=c(2.5,4,0.1,1),my_inset=c(0,-0.2))
#my_ylim = c(logodds100(10.5),logodds100(40))
text(1992, logodds100(4), "B. Natives",cex = .8)
if (f2$ltr_all$coefficients["year"]>0){plus <- "+"
}  else {plus <- ""}
my_delta <- paste("\n",plus,round(f2$ltr_all$coefficients["year"],1),"% /year",sep="")
text(f2$des$mystart[f2$des$country %in% "Adj. Mean"],
     logodds100(-0.02+f2$des$first[f2$des$country %in% "Adj. Mean"]), 
     bquote (Delta : .(my_delta)),cex = .6,pos=1)
dev.off()

## Figure S5.2 Migrants and female workplace relative isolation ####

# Import gdr data
seg_gdr <- readRDS("Data/rds/exp_seg_gdr.rds")

seg_gdr <- seg_gdr[(seg_gdr$country=="Denmark")+(seg_gdr$year==1994)<2,]
seg_gdr1 <- seg_gdr[seg_gdr$female %in% 1,]
seg_gdr0 <- seg_gdr[seg_gdr$female %in% 0,]
seg_gdr1$ofemale <- oddsratio(seg_gdr1$mean_xfemale,seg_gdr0$mean_xfemale)
seg_gdr1$lofemale <- log(seg_gdr1$ofemale)

#Figure
pdf(file="Figures/Figure_S5_2.pdf",width=7,height=10,family="serif")
par(mfrow=c(2,1))
mig_eo <- mygraph_f(seg_mig1$lomigrant,"",
                    "Relative exposure (odds-ratio - log scale)","",
                    ratio="logodds",legpos="top",
                    myyaxt="n",
                    my_mar=c(0.1,4,3,1),my_inset=c(0,-0.1),myxaxt="n",my_legend = "n",my_xlim=c(1990,2019))
text(1991, log(5), "A. Migrant",cex = .8)
if (mig_eo$ltr_all$coefficients["year"]>0){plus <- "+"
}  else {plus <- ""}
my_delta <- paste("\n",plus,round(mig_eo$ltr_all$coefficients["year"],1),"% /year",sep="")
text(mig_eo$des$mystart[mig_eo$des$country %in% "Adj. Mean"],
     -0.01+mig_eo$des$first[mig_eo$des$country %in% "Adj. Mean"], 
     bquote (Delta : .(my_delta)),cex = .6,pos=1)
legend("top", 
       inset=c(0,-0.15),
       legend=c("Adj. Mean",myctry),pch =c(NA,mypch), 
       col=c("black",mycol),
       pt.bg=c("black",mycol),
       cex=0.8,ncol=6,lty=c(1,mylty),lwd = c(3,rep(1,length(myctry))),
       text.col = c("black",mycol),bg=NA,box.col="gray",box.lty=0)
gdr_eo <- mygraph_f(seg_gdr1$lofemale,"",
                    "Relative exposure (odds-ratio - log scale)","",
                    ratio="logodds",legpos="top",
                    myyaxt="n",my_legend = "n",
                    my_mar=c(2.5,4,0.1,1),my_inset=c(0,-0.2))
#my_ylim = c(logodds100(10.5),logodds100(40))

if (gdr_eo$ltr_all$coefficients["year"]>0){plus <- "+"
}  else {plus <- ""}
my_delta <- paste("\n",plus,round(gdr_eo$ltr_all$coefficients["year"],1),"% /year",sep="")

text(1991, log(4), "B. Female",cex = .8)

text(gdr_eo$des$mystart[gdr_eo$des$country %in% "Adj. Mean"],
     -0.01+gdr_eo$des$first[gdr_eo$des$country %in% "Adj. Mean"], 
     bquote (Delta : .(my_delta)),cex = .6,pos=1)
dev.off()

## Figure S5.3. Age isolation ####

# Data
age <- readRDS("Data/rds/exp_seg_age.rds")

colnames(age)[which(names(age) %in% c("mean_xage1","mean_xage2","mean_xage3","mean_xage4"))] <- c("xage1","xage2","xage3","xage4")
age$xage_t <- age$xage1+age$xage2+age$xage3+age$xage4
age$xage1 <- age$xage1/age$xage_t
age$xage2 <- age$xage2/age$xage_t
age$xage3 <- age$xage3/age$xage_t
age$xage4 <- age$xage4/age$xage_t
age$xage_t <- age$xage1+age$xage2+age$xage3+age$xage4  

age_p <- age[age$age %in% NA,c("year","country","xage1","xage2","xage3","xage4")]
colnames(age_p)[which(names(age_p) %in% c("xage1","xage2","xage3","xage4"))] <- c("age1","age2","age3","age4")
age_p$t <- age_p$age1+age_p$age2+age_p$age3+age_p$age4
age <- merge(age,age_p,by=c("year","country"),all.x=TRUE)

age_w1 <- age[age$age %in% 1,c("year","country","sum_w_xage1")]
age_w2 <- age[age$age %in% 2,c("year","country","sum_w_xage1")]
age_w3 <- age[age$age %in% 3,c("year","country","sum_w_xage1")]
age_w4 <- age[age$age %in% 4,c("year","country","sum_w_xage1")]
age_wall <- age[age$age %in% NA,c("year","country","sum_w_xage1")]

age_w <- merge(age_w1,age_w2,by=c("year","country"),all=TRUE)
colnames(age_w) <- c("year","country","wgtage1","wgtage2")
age_w <- merge(age_w,age_w3,by=c("year","country"),all=TRUE)
colnames(age_w) <- c("year","country","wgtage1","wgtage2","wgtage3")
age_w <- merge(age_w,age_w4,by=c("year","country"),all=TRUE)
colnames(age_w) <- c("year","country","wgtage1","wgtage2","wgtage3","wgtage4")
age_w <- merge(age_w,age_wall,by=c("year","country"),all=TRUE)
colnames(age_w) <- c("year","country","wgtage1","wgtage2","wgtage3","wgtage4","swgt")  

age_w$wgtage1 <- age_w$wgtage1/age_w$swgt
age_w$wgtage2 <- age_w$wgtage2/age_w$swgt
age_w$wgtage3 <- age_w$wgtage3/age_w$swgt
age_w$wgtage4 <- age_w$wgtage4/age_w$swgt

age <- merge(age,age_w,by=c("year","country"),all.x=TRUE)

str(age)
age$xage1_n <- ifelse(age$age==1,(age$age1-age$wgtage1*age$xage1)/(1-age$wgtage1),
                      ifelse(age$age==2,(age$age1-age$wgtage2*age$xage1)/(1-age$wgtage2),
                             ifelse(age$age==3,(age$age1-age$wgtage3*age$xage1)/(1-age$wgtage3),
                                    ifelse(age$age==4,(age$age1-age$wgtage4*age$xage1)/(1-age$wgtage4),NA))))

age$xage2_n <- ifelse(age$age==1,(age$age2-age$wgtage1*age$xage2)/(1-age$wgtage1),
                      ifelse(age$age==2,(age$age2-age$wgtage2*age$xage2)/(1-age$wgtage2),
                             ifelse(age$age==3,(age$age2-age$wgtage3*age$xage2)/(1-age$wgtage3),
                                    ifelse(age$age==4,(age$age2-age$wgtage4*age$xage2)/(1-age$wgtage4),NA))))


age$xage3_n <- ifelse(age$age==1,(age$age3-age$wgtage1*age$xage3)/(1-age$wgtage1),
                      ifelse(age$age==2,(age$age3-age$wgtage2*age$xage3)/(1-age$wgtage2),
                             ifelse(age$age==3,(age$age3-age$wgtage3*age$xage3)/(1-age$wgtage3),
                                    ifelse(age$age==4,(age$age3-age$wgtage4*age$xage3)/(1-age$wgtage4),NA))))                                  

age$xage4_n <- ifelse(age$age==1,(age$age4-age$wgtage1*age$xage4)/(1-age$wgtage1),
                      ifelse(age$age==2,(age$age4-age$wgtage2*age$xage4)/(1-age$wgtage2),
                             ifelse(age$age==3,(age$age4-age$wgtage3*age$xage4)/(1-age$wgtage3),
                                    ifelse(age$age==4,(age$age4-age$wgtage4*age$xage4)/(1-age$wgtage4),NA))))                                  




age$odage1 <- ((age$xage1/(1-age$xage1))/((age$xage1_n/(1-age$xage1_n))))
age$odage2 <- ((age$xage2/(1-age$xage2))/((age$xage2_n/(1-age$xage2_n))))
age$odage3 <- ((age$xage3/(1-age$xage3))/((age$xage3_n/(1-age$xage3_n))))
age$odage4 <- ((age$xage4/(1-age$xage4))/((age$xage4_n/(1-age$xage4_n))))


age$lodage1 <- log(age$odage1)
age$lodage2 <- log(age$odage2)
age$lodage3 <- log(age$odage3)
age$lodage4 <- log(age$odage4)


seg_age4 <- age[age$age==4,]
seg_age1 <- age[age$age==1,]
 

# Figure
pdf(file="Figures/Figure_S5_3.pdf",width=7,height=10,family="serif")
par(mfrow=c(2,1))
f1 <- mygraph_f(seg_age4$lodage4,"",
                "Relative exposure (odds-ratio - log scale)","",
                ratio="logodds",legpos="top",
                myyaxt="n",
                my_mar=c(0.1,4,3,1),my_inset=c(0,-0.1),myxaxt="n",my_legend = "n")
text(1991, log(2.5), "A. Age>55",cex = .8)
if (f1$ltr_all$coefficients["year"]>0){plus <- "+"
}  else {plus <- ""}
my_delta <- paste("\n",plus,round(f1$ltr_all$coefficients["year"],1),"% /year",sep="")
text(f1$des$mystart[f1$des$country %in% "Adj. Mean"],
     -0.05+f1$des$first[f1$des$country %in% "Adj. Mean"], 
     bquote (Delta : .(my_delta)),cex = .6,pos=1)
legend("top", 
       inset=c(0,-0.15),
       legend=c("Adj. Mean",myctry[c(1:13)]),
       pch =c(NA,mypch[c(1:13)]), 
       col=c("black",mycol[c(1:13)]),
       pt.bg=c("black",mycol[c(1:13)]),
       lty=c(1,mylty[c(1:13)]),
       lwd = c(3,rep(1,length(myctry[c(1:13)]))),
       text.col = c("black",mycol[c(1:13)]),
       cex=0.8,ncol=6,
       bg=NA,box.col="gray",box.lty=0)
f2 <- mygraph_f(seg_age1$lodage1,"",
                "Relative exposure (odds-ratio - log scale)","",
                ratio="logodds",legpos="top",
                myyaxt="n",my_legend = "n",
                my_mar=c(2.5,4,0.1,1),my_inset=c(0,-0.2))
#my_ylim = c(logodds100(10.5),logodds100(40))
text(1991, log(2.5), "B. Age<31",cex = .8)
if (f2$ltr_all$coefficients["year"]>0){plus <- "+"
}  else {plus <- ""}
my_delta <- paste("\n",plus,round(f2$ltr_all$coefficients["year"],1),"% /year",sep="")
text(f2$des$mystart[f2$des$country %in% "Adj. Mean"],
     -0.01+f2$des$first[f2$des$country %in% "Adj. Mean"], 
     bquote (Delta : .(my_delta)),cex = .6,pos=1)
dev.off()

## Figure S5.4. Class isolation ####

# Data occupation
occ <- readRDS("Data/rds/exp_seg_occ.rds")

colnames(occ)[which(names(occ) %in% c("mean_xocc1","mean_xocc2","mean_xocc3"))] <- c("xocc1","xocc2","xocc3")

occ_p <- occ[occ$occ %in% NA,c("year","country","xocc1","xocc2","xocc3")]
colnames(occ_p)[which(names(occ_p) %in% c("xocc1","xocc2","xocc3"))] <- c("occ1","occ2","occ3")
occ_p$t <- occ_p$occ1+occ_p$occ2+occ_p$occ3
occ <- merge(occ,occ_p,by=c("year","country"),all.x=TRUE)

occ_w1 <- occ[occ$occ %in% 1,c("year","country","sum_w_xocc1")]
occ_w2 <- occ[occ$occ %in% 2,c("sum_w_xocc1")]
occ_w3 <- occ[occ$occ %in% 3,c("sum_w_xocc1")]
occ_wall <- occ[is.na(occ$occ),c("sum_w_xocc1")]

rownames(occ_w1) <- NULL
occ_w <- cbind(occ_w1,occ_w2,occ_w3,occ_wall)  
colnames(occ_w) <- c("year","country","wgtocc1","wgtocc2","wgtocc3","swgt")  

occ_w$wgtocc1 <- occ_w$wgtocc1/occ_w$swgt
occ_w$wgtocc2 <- occ_w$wgtocc2/occ_w$swgt
occ_w$wgtocc3 <- occ_w$wgtocc3/occ_w$swgt

occ <- merge(occ,occ_w,by=c("year","country"),all.x=TRUE)

occ$xocc1_n <- ifelse(occ$occ==1,(occ$occ1-occ$wgtocc1*occ$xocc1)/(1-occ$wgtocc1),
                      ifelse(occ$occ==2,(occ$occ1-occ$wgtocc2*occ$xocc1)/(1-occ$wgtocc2),
                             ifelse(occ$occ==3,(occ$occ1-occ$wgtocc3*occ$xocc1)/(1-occ$wgtocc3),NA)))

occ$xocc2_n <- ifelse(occ$occ==1,(occ$occ2-occ$wgtocc1*occ$xocc2)/(1-occ$wgtocc1),
                      ifelse(occ$occ==2,(occ$occ2-occ$wgtocc2*occ$xocc2)/(1-occ$wgtocc2),
                             ifelse(occ$occ==3,(occ$occ2-occ$wgtocc3*occ$xocc2)/(1-occ$wgtocc3),NA)))

occ$xocc3_n <- ifelse(occ$occ==1,(occ$occ3-occ$wgtocc1*occ$xocc3)/(1-occ$wgtocc1),
                      ifelse(occ$occ==2,(occ$occ3-occ$wgtocc2*occ$xocc3)/(1-occ$wgtocc2),
                             ifelse(occ$occ==3,(occ$occ3-occ$wgtocc3*occ$xocc3)/(1-occ$wgtocc3),NA)))


occ$odocc1 <- ((occ$xocc1/(1-occ$xocc1))/((occ$xocc1_n/(1-occ$xocc1_n))))
occ$odocc2 <- ((occ$xocc2/(1-occ$xocc2))/((occ$xocc2_n/(1-occ$xocc2_n))))
occ$odocc3 <- ((occ$xocc3/(1-occ$xocc3))/((occ$xocc3_n/(1-occ$xocc3_n))))


occ$lodocc1 <- log(occ$odocc1)
occ$lodocc2 <- log(occ$odocc2)
occ$lodocc3 <- log(occ$odocc3)

occ$docc1 <- (occ$xocc1-occ$xocc1_n)
occ$docc2 <- (occ$xocc2-occ$xocc2_n)
occ$docc3 <- (occ$xocc3-occ$xocc3_n)


seg_occna <- occ[is.na(occ$occ),]

seg_occ1 <- occ[occ$occ %in% 1
                & ((occ$country=="Denmark")+(occ$year==1994 | occ$year>2009))<2
                & ((occ$country=="Germany")+(occ$year>=2011))<2
                & ((occ$country=="Korea")+(occ$year>=2008))<2
                & ((occ$country=="Korea")+(occ$year<=1992))<2
                ,]
seg_occ3 <- occ[occ$occ %in% 3
                & ((occ$country=="Denmark")+(occ$year==1994 | occ$year>2009))<2
                & ((occ$country=="Germany")+(occ$year>=2011))<2
                & ((occ$country=="Korea")+(occ$year>=2008))<2
                & ((occ$country=="Korea")+(occ$year<=1992))<2
                ,]

# Figure
pdf(file="Figures/Figure_S5_4.pdf",width=7,height=10,family="serif")
par(mfrow=c(2,1))
f1 <- mygraph_f(seg_occ1$lodocc1,"",
                "Relative exposure (odds-ratio - log scale)","",
                ratio="logodds",legpos="top",
                myyaxt="n",
                my_mar=c(0.1,4,3,1),my_inset=c(0,-0.1),myxaxt="n",my_legend = "n")
text(1995, log(9), "A. Managers and professionals",cex = .8)
if (f1$ltr_all$coefficients["year"]>0){plus <- "+"
}  else {plus <- ""}
my_delta <- paste("\n",plus,round(f1$ltr_all$coefficients["year"],1),"% /year",sep="")
text(f1$des$mystart[f1$des$country %in% "Adj. Mean"],
     -0.05+f1$des$first[f1$des$country %in% "Adj. Mean"], 
     bquote (Delta : .(my_delta)),cex = .6,pos=1)
field <- c(NA,2:6,8,NA,9:13)
legend("top", 
       inset=c(0-0.15),
       legend=c("Adj. Mean",myctry[field]),
       pch =c(NA,mypch[field]), 
       col=c("black",mycol[field]),
       pt.bg=c("black",mycol[field]),
       lty=c(1,mylty[field]),
       lwd = c(3,rep(1,length(myctry[field]))),
       text.col = c("black",mycol[field]),
       cex=0.8,ncol=6,
       bg=NA,box.col="gray",box.lty=0)
f2 <- mygraph_f(seg_occ3$lodocc3,"",
                "Relative exposure (odds-ratio - log scale)","",
                ratio="logodds",legpos="top",
                myyaxt="n",my_legend = "n",
                my_mar=c(2.5,4,0.1,1),my_inset=c(0,-0.2))
text(1993, log(8), "B. Working-class employees",cex = .8)
if (f2$ltr_all$coefficients["year"]>0){plus <- "+"
}  else {plus <- ""}
my_delta <- paste("\n",plus,round(f2$ltr_all$coefficients["year"],1),"% /year",sep="")
text(f2$des$mystart[f2$des$country %in% "Adj. Mean"],
     -0.01+f2$des$first[f2$des$country %in% "Adj. Mean"], 
     bquote (Delta : .(my_delta)),cex = .6,pos=1)
dev.off()


## Figure S9.1 ####

fr_ev <- readRDS("Data/rds/exp_fr_ev.rds")
dd_out <- fr_ev[fr_ev$event=="out2",]

# Outsourcing 
plot_out <- ggplot(dd_out, aes(x=time, y=b,color=ctrl)) +        # ggplot2 plot with confidence intervals
  theme_bw()+
  theme(panel.grid.major = element_line(color = "white"))+
  geom_hline(yintercept=0,color="black")+
  ylab("Effect of 10 percent workforce outsourcing on top 10% isolation (in percentage point)") +
  geom_vline(xintercept=-0.5,color="black",linetype="dashed",size=1)+
  # geom_vline(xintercept=2011.5,color="gray5",linetype="dashed")+
  geom_point(aes(shape=ctrl),size=3) +
  geom_line() +
  geom_errorbar(aes(ymin = l, ymax = u),width=0.5)+
  scale_x_continuous(breaks = round(seq(min(floor(dd_out$time)), max(dd_out$time), by = 1),1))+
  # annotate("text", x=2010, y=-2.0, label= "Offshoring event
  #         between 2009 & 2011",cex=3) +
  # annotate("text", x=2006.5, y=5, label= "Pre-event difference",col="red") +
  # annotate("text", x=2013.5, y=-2, label= "Post-event difference",col="blue") +
  theme(legend.position="bottom")+ 
  theme(legend.title = element_blank())+
  ylab("Top 10% isolation (in % point)") +
  ggtitle("Outsourcing event")+
  scale_color_manual(values=c("darkblue", "darkorange"))+
  theme(legend.position="none") +
  # theme(axis.title.y = element_blank()) +
  theme(axis.title.x = element_blank())
dev.off()
plot_out


#-------------------------------------------------#
#### Layoff ####

dd_lay <- fr_ev[fr_ev$event=="lay2",]
# Plot 
plot_lay <- ggplot(dd_lay, aes(x=time, y=b,color=ctrl)) +        # ggplot2 plot with confidence intervals
  theme_bw()+
  theme(panel.grid.major = element_line(color = "white"))+
  geom_hline(yintercept=0,color="black")+
  geom_vline(xintercept=-0.5,color="black",linetype="dashed",size=1)+
  # geom_vline(xintercept=2011.5,color="gray5",linetype="dashed")+
  geom_point(aes(shape=ctrl),size=3) +
  geom_line() +
  geom_errorbar(aes(ymin = l, ymax = u),width=0.5)+
  scale_x_continuous(breaks = round(seq(min(floor(dd_lay$time)), max(dd_lay$time), by = 1),1))+
  # annotate("text", x=2010, y=-2.0, label= "Offshoring event
  #         between 2009 & 2011",cex=3) +
  # annotate("text", x=2006.5, y=5, label= "Pre-event difference",col="red") +
  # annotate("text", x=2013.5, y=-2, label= "Post-event difference",col="blue") +
  theme(legend.position="bottom")+
  theme(legend.title = element_blank())+
  theme(legend.position="none") +
  scale_color_manual(values=c("darkblue", "darkorange"))+
  ylab("Top 10% isolation (in % point)") +
  # theme(axis.title.y = element_blank()) +
  # theme(axis.title.x = element_blank()) +
  
  ggtitle("Layoff event") 
dev.off()
plot_lay


#-------------------------------------------------#
#### Offshoring ####
dd_off <- fr_ev[fr_ev$event=="off2",]

# Plot 
plot_off <- ggplot(dd_off, aes(x=time, y=b,color=ctrl)) +        # ggplot2 plot with confidence intervals
  theme_bw()+
  theme(panel.grid.major = element_line(color = "white"))+
  geom_hline(yintercept=0,color="black")+
  # ylab("Effect of 10 percent workforce offshoring on top 10% isolation (in percentage point)") +
  geom_vline(xintercept=2008.5,color="black",linetype="dashed",size=1)+
  geom_vline(xintercept=2011.5,color="gray5",linetype="dashed")+
  geom_point(aes(shape=ctrl),size=3) +
  geom_line() +
  geom_errorbar(aes(ymin = l, ymax = u),width=0.5)+
  scale_x_continuous(breaks = round(seq(min(floor(dd_off$time)), max(dd_off$time), by = 1),1))+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  annotate("text", x=2010, y=-0.4, label= "Event between
  2009 & 2011",cex=3) +
  # annotate("text", x=2006.5, y=5, label= "Pre-event difference",col="red") +
  # annotate("text", x=2013.5, y=-2, label= "Post-event difference",col="blue") +
  theme(legend.position="bottom")+ 
  theme(legend.title = element_blank())+
  ylab("Top 10% isolation (in % point)") +
  ggtitle("Offshoring event")+
  theme(legend.position="none") +
  scale_color_manual(values=c("darkblue", "darkorange"))+
  theme(axis.title.y = element_blank()) +
  theme(axis.title.x = element_blank())
dev.off()
plot_off

#### Subcontracting ####
dd_sub <- fr_ev[fr_ev$event=="sub2",]

# Plot 
plot_sub <- ggplot(dd_sub, aes(x=time, y=b,color=ctrl)) +        # ggplot2 plot with confidence intervals
  theme_bw()+
  theme(panel.grid.major = element_line(color = "white"))+
  geom_hline(yintercept=0,color="black")+
  # ylab("Effect of 10 percent workforce offshoring on top 10% isolation (in percentage point)") +
  geom_vline(xintercept=2005.5,color="black",linetype="dashed",size=1)+
  geom_vline(xintercept=2010.5,color="gray5",linetype="dashed")+
  geom_point(aes(shape=ctrl),size=3) +
  geom_line() +
  geom_errorbar(aes(ymin = l, ymax = u),width=0.5)+
  scale_x_continuous(breaks = round(seq(min(floor(dd_sub$time)), max(dd_sub$time), by = 1),1))+
  annotate("text", x=2008, y=-0.5, label= "Event between 2005 & 2011",cex=3) +
  # annotate("text", x=2006.5, y=5, label= "Pre-event difference",col="red") +
  # annotate("text", x=2013.5, y=-2, label= "Post-event difference",col="blue") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  theme(legend.position="bottom")+ 
  theme(legend.title = element_blank())+
  ylab("Top 10% isolation (in % point)") +
  ggtitle("Subcontracting event")+
  # theme(legend.position="none") +
  scale_color_manual(values=c("darkblue", "darkorange"))+
  theme(axis.title.y = element_blank()) +
  theme(axis.title.x = element_blank())
dev.off()
plot_sub



grid.arrange(plot_out, plot_off, plot_lay, plot_sub, ncol=2)

pdf("Figures/Figure_S91_France_event.pdf")
grid.arrange(plot_out, plot_off, plot_lay, plot_sub, ncol=2)
dev.off()



#---------------------------------#
#---------------------------------#
#---------------------------------#
# TABLES ####
#---------------------------------#
#---------------------------------#
#---------------------------------#


#---------------------------------#
## Table 1. ELEMENTS FOR TABLE 1 ####
#---------------------------------#

all_e <- readRDS("Data/rds/exp_all_e.rds")
all_f <- readRDS("Data/rds/exp_all_f.rds")

all_e$nb_units <- all_e$mean_invcountuny*all_e$sum_w_countuny
all_f$nb_units <- all_f$mean_invcountuny*all_f$sum_w_countuny




minyear <- tapply(all_e$year,all_e$country,min)
maxyear <- tapply(all_e$year,all_e$country,max)
maxyear-minyear+1

minwage <- tapply(all_e$min_wage,list(all_e$year,all_e$country),min,na.rm=T)
nbwork <- tapply(all_e$nb_obs[is.na(all_e$fwage) & all_e$year>=2012],list(all_e$year[is.na(all_e$fwage) & all_e$year>=2012],all_e$country[is.na(all_e$fwage) & all_e$year>=2012]),min,na.rm=T)
nbunits <- tapply(all_e$nb_units[is.na(all_e$fwage) & all_e$year>=2012],list(all_e$year[is.na(all_e$fwage) & all_e$year>=2012],all_e$country[is.na(all_e$fwage) & all_e$year>=2012]),min,na.rm=T)
nbunits_f <- tapply(all_f$nb_units[is.na(all_f$fwage) & all_f$year>=2012],list(all_f$year[is.na(all_f$fwage) & all_f$year>=2012],all_f$country[is.na(all_f$fwage) & all_f$year>=2012]),min,na.rm=T)


all_e2 <- all_e
all_e2$maxyear <- ave(all_e$year,all_e$country,FUN=max)
all_e2 <- all_e2[all_e2$year ==all_e2$maxyear & is.na(all_e2$fwage), ]

all_f2 <- all_f
all_f2$maxyear <- ave(all_f$year,all_f$country,FUN=max)
all_f2 <- all_f2[all_f2$year ==all_f2$maxyear & is.na(all_f2$fwage), ]


minwage2 <- tapply(all_e2$min_wage,all_e2$country,mean,na.rm=T)
nbwork2 <- tapply(all_e2$nb_obs,all_e2$country,min,na.rm=T)
nbunits2 <- tapply(all_e2$nb_units,all_e2$country,min,na.rm=T)
nbunits_f2 <- tapply(all_f2$nb_units,all_f2$country,min,na.rm=T)

table_1 <- cbind(minyear,maxyear,minwage2,round(nbwork2,0),round(nbunits2,0))
table_1 <- merge(table_1,round(nbunits_f2,0),by='row.names',all.x=T)
colnames(table_1) <- c("country","start","end","end_minwage","end_nb_workers","end_nb_est","end_nb_firm")
table_1

# Czechia and Netherlands data production did not allow to estimate mins
write.csv(table_1,"Tables/Table 1.csv",row.names = F,na=".")


#---------------------------------#---------------------------------#
## Table 2. Main Trends ####
#---------------------------------#---------------------------------#

#Preliminary figures for producing statistics

  all_e <- readRDS("Data/rds/exp_all_e.rds")
  all_e5 <- all_e[all_e$fwage %in% 5,]
  all_e45 <- all_e[all_e$fwage %in% 45,]
  all_e1 <- all_e[all_e$fwage %in% 1,]
  
  d_all <- readRDS("Data/rds/exp_d_all.rds")
  d0 <- d_all[d_all$d %in% 0,]
  
  all_e$var_within <- all_e$sd_lnwage_within**2
  all_e$var_total <- all_e$sd_lnwage**2
  all_e$var_between <- all_e$var_total-all_e$var_within
  all_e$pvar_within <- all_e$var_within/all_e$var_total
  all_e$pvar_between <- all_e$var_between/all_e$var_total
  all_e$nb_units <- all_e$mean_invcountuny*all_e$sum_w_countuny
  all_evarlwage <- all_e[all_e$fwage %in% NA,c("year","country","var_within","var_total","var_between","pvar_within","pvar_between")]
  
  top1share <- 0.01*all_e[all_e$fwage %in% 5,"mean_wage"]/all_e[all_e$fwage %in% NA,"mean_wage"]
  top10share <- 0.1*(0.9*all_e[all_e$fwage %in% 4,"mean_wage"]+0.1*all_e[all_e$fwage %in% 5,"mean_wage"])/all_e[all_e$fwage %in% NA,"mean_wage"]
  all_eshare <- data.frame(top1share,top10share,all_e[all_e$fwage %in% NA,c("mean_wage","country","year")])
  
  
  g_e55 <- mygraph2(all_e5$mean_xF9910,"Figure 1A. Top 1% workplace isolation",
                    "Exposure (% - log-odds scale)","",ratio="OR",legpos="bottomright",myyaxt="n")
  axis(2,logodds100(c(5,11,13,15,17,19,21,22,23)),label = paste(c(5,11,13,15,17,19,21,22,23),"%"), las=2,cex.axis=0.8)
  
  
  g_e4545 <- mygraph2(all_e45$mean_xF9010,"Figure 1B. Top 10% workplace isolation",
                      "Exposure (% - log-odds scale)","",ratio="OR",legpos = "topleft",myyaxt="n")
  axis(2,logodds100(c(21:24,26:29,31:34,36:39,41:44)),label = paste(c(21:24,26:29,31:34,36:39,41:44),"%"), las=2,cex.axis=0.8)
  
  #Figure S2A
  g_e52 <- mygraph2(all_e5$mean_xF2575,"Figure S1A. Top 1%  workplace exposure to mid quartiles","Exposure (% - log-odds scale)","",ratio="OR",legpos="bottomleft",myyaxt="n")
  
  
  #Figure S2B
  g_e452 <- mygraph2(all_e45$mean_xF2575,"Figure S1B. Top 10%  workplace exposure to mid quartiles","Exposure (% - log-odds scale)","",ratio="OR",legpos="bottomleft",myyaxt="n")
  
  
  g_e51 <- mygraph2(all_e5$mean_xF0025,"Figure 2A. Top 1% workplace exposure to bottom 25%",
                    "Exposure (% - log-odds scale)","",ratio="OR",legpos = "bottomleft",myyaxt="n")
  #matlines(g_e51$g$year,rep(logodds100(9),29),lwd = 1,col="gray")
  
  g_e451 <- mygraph2(all_e45$mean_xF0025,"Figure 2B. Top 10% workplace exposure to bottom 25%",
                     "Exposure (% - log-odds scale)","",ratio="OR",legpos = "topright",myyaxt="n")
  axis(2,logodds100(c(11,13,15)),label = paste(c(11,13,15),"%"), las=2,cex.axis=0.8)
  
  
  g_e11 <- mygraph2(all_e1$mean_xF0025,"Figure 3. Bottom 25%  workplace isolation","Exposure (% - log-odds scale)","",ratio="OR",legpos="topleft",myyaxt="n")
  #axis(2,logodds100(c(36:39,41:44,46:49,51:54,56:59,61:64)),label = paste(c(36:39,41:44,46:49,51:54,56:59,61:64),"%"), las=2,cex.axis=0.8)
  
  g_d0 <- mygraph2(d0$mean_xd0,"Bottom 10%  isolation","Exposure (% - log-odds scale)","",ratio="OR",legpos="topright",myyaxt="n")
  
  beshlv <- mygraph2(all_evarlwage$pvar_between,"Between share of log wage variance",
                     "Share (% - log-odds scale)","",ratio="OR",legpos = "topleft",myyaxt="n")
  
  sh1 <- mygraph2(all_eshare$top1share,"Top 1% wage share","Share of wage bill (% - log-odds scale)","",ratio="OR",legpos = "topleft",myyaxt="n")
  axis(2,logodds100(c(11)),label = paste(c(11),"%"), las=2,cex.axis=0.8)
  
  sh10 <- mygraph2(all_eshare$top10share,"Top 10% wage share","Share of wage bill (% - log-odds scale)","",ratio="OR",legpos = "topleft",myyaxt="n")
  
  dev.off()

#END YEAR LEVEL
tab0 <- rbind.fill(data.frame(t(data.frame(g_e55$des[,c("last","gr_rate","myend")])),c("last","gr_rate","kend"),c("top1isol","top1isol","top1isol")),
                   data.frame(t(g_e4545$des[,c("last","gr_rate")]),c("last","gr_rate"),c("top10isol","top10isol")),
                   data.frame(t(g_e52$des[,c("last","gr_rate")]),c("last","gr_rate"),c("top1_midq","top1_midq")),
                   data.frame(t(g_e452$des[,c("last","gr_rate")]),c("last","gr_rate"),c("top10_midq","top10_midq")),
                   data.frame(t(g_e51$des[,c("last","gr_rate")]),c("last","gr_rate"),c("top1_bottomq","top1_bottomq")),
                   data.frame(t(g_e451$des[,c("last","gr_rate")]),c("last","gr_rate"),c("top10_bottomq","top10_bottomq")),
                   data.frame(t(g_e11$des[,c("last","gr_rate")]),c("last","gr_rate"),c("bottom25isol","bottom25isol")),
                   data.frame(t(g_d0$des[,c("last","gr_rate")]),c("last","gr_rate"),c("bottom10isol","bottom10isol")),
                   data.frame(t(beshlv$des[,c("last","gr_rate")]),c("last","gr_rate"),c("betw_sh","betw_sh")),
                   data.frame(t(sh1$des[,c("last","gr_rate")]),c("last","gr_rate"),c("sh1","sh1")),
                   data.frame(t(sh10$des[,c("last","gr_rate")]),c("last","gr_rate"),c("sh10","sh10")))

str(tab0)
tab0$type <- gsub("NA","",paste(as.character(tab0[,15]),as.character(tab0[,17]),sep=""))
tab0$nature <- gsub("NA","",paste(
  as.character(tab0[,16]),
  as.character(tab0[,18]),
  as.character(tab0[,19]),
  as.character(tab0[,20]),
  as.character(tab0[,21]),
  as.character(tab0[,22]),
  as.character(tab0[,23]),
  as.character(tab0[,24]),
  as.character(tab0[,25]),
  as.character(tab0[,26]),
  as.character(tab0[,27]),
  sep=""))

str(tab0)
tab1 <- tab0[order(tab0$type) ,c("type","nature","Canada","Denmark","Norway","Sweden","France","Netherlands","Germany","Spain",
                                 "Czechia","Hungary","Korea","Japan","Adj..Mean")]

table(tab0$type)
table_2A <- t(tab1[tab1$type %in% c("kend","last"),])
write.csv(table_2A,"Tables/Table 2A.csv",row.names = F,na=".")


#TRENDS
m_e55<- felm(I(100*logodds(mean_xF9910))~country:year|country|0|year,
             data=all_e5[(all_e5$country %in% "France" & all_e5$year==1993)==F
                         & (all_e5$country %in% "Denmark" & all_e5$year==1994)==F
                         ,])
screenreg(m_e55)

m_e4545<- felm(I(100*logodds(mean_xF9010))~country:year|country|0|year,
               data=all_e45[(all_e45$country %in% "France" & all_e45$year==1993)==F 
                            & (all_e45$country %in% "Denmark" & all_e45$year==1994)==F
                            ,])
screenreg(m_e4545)

m_e52<- felm(I(100*logodds(mean_xF2575))~country:year|country|0|year,
             data=all_e5[(all_e5$country %in% "France" & all_e5$year==1993)==F  
                         & (all_e5$country %in% "Denmark" & all_e5$year==1994)==F
                         ,])
screenreg(m_e52)

m_e452<- felm(I(100*logodds(mean_xF2575))~country:year|country|0|year,
              data=all_e45[(all_e45$country %in% "France" & all_e45$year==1993)==F  
                           & (all_e45$country %in% "Denmark" & all_e45$year==1994)==F
                           ,])
screenreg(m_e452)


m_e51 <- felm(I(100*logodds(mean_xF0025))~country:year|country|0|year,
              data=all_e5[(all_e5$country %in% "France" & all_e5$year==1993)==F  
                          & (all_e5$country %in% "Denmark" & all_e5$year==1994)==F
                          ,])
screenreg(m_e51)

m_e451 <- felm(I(100*logodds(mean_xF0025))~country:year|country|0|year,
               data=all_e45[
                              (all_e45$country %in% "France" & all_e45$year==1993)==F
                            & (all_e45$country %in% "Denmark" & all_e45$year==1994)==F
                            ,])
screenreg(m_e451)


m_e11 <- felm(I(100*logodds(mean_xF0025))~country:year|country|0|year,
              data=all_e1[(all_e1$country %in% "France" & all_e1$year==1993)==F
                          & (all_e1$country %in% "Denmark" & all_e1$year==1994)==F
                          ,])
screenreg(m_e11)


m_d0 <- felm(I(100*logodds(mean_xd0))~country:year|country|0|year,
             data=d0[(d0$country %in% "France" & d0$year==1993)==F
                     & (d0$country %in% "Denmark" & d0$year==1994)==F
                     ,])
screenreg(m_d0)

m_beshlv <- felm(I(100*logodds(pvar_between))~country:year|country|0|year,
                 data=all_evarlwage[(all_evarlwage$country %in% "France" & all_evarlwage$year==1993)==F
                                    & (all_evarlwage$country %in% "Denmark" & all_evarlwage$year==1994)==F
                                    ,])
screenreg(m_beshlv)


m_sh1 <- felm(I(100*logodds(top1share))~country:year|country|0|year,
              data=all_eshare[(all_eshare$country %in% "France" & all_eshare$year==1993)==F
                              & (all_eshare$country %in% "Denmark" & all_eshare$year==1994)==F
                              ,])
screenreg(m_sh1)

m_sh10 <- felm(I(100*logodds(top10share))~country:year|country|0|year,
               data=all_eshare[(all_eshare$country %in% "France" & all_eshare$year==1993)==F
                               & (all_eshare$country %in% "Denmark" & all_eshare$year==1994)==F
                               ,])
screenreg(m_sh10)


m_e55_all <-  felm(I(100*logodds(mean_xF9910))~year|country|0|year,
                   data=all_e5[(all_e5$country %in% "France" & all_e5$year==1993)==F
                               & (all_e5$country %in% "Denmark" & all_e5$year==1994)==F
                               ,])
screenreg(m_e55)

m_e4545_all <-  felm(I(100*logodds(mean_xF9010))~year|country|0|year,
                     data=all_e45[(all_e45$country %in% "France" & all_e45$year==1993)==F
                                  & (all_e45$country %in% "Denmark" & all_e45$year==1994)==F
                                  ,])
screenreg(m_e4545)

m_e52_all <-  felm(I(100*logodds(mean_xF2575))~year|country|0|year,
                   data=all_e5[
                                 (all_e5$country %in% "France" & all_e5$year==1993)==F
                               & (all_e5$country %in% "Denmark" & all_e5$year==1994)==F
                               ,])
screenreg(m_e52)

m_e452_all <-  felm(I(100*logodds(mean_xF2575))~year|country|0|year,
                    data=all_e45[
                                   (all_e45$country %in% "France" & all_e45$year==1993)==F
                                 & (all_e45$country %in% "Denmark" & all_e45$year==1994)==F
                                 ,])
screenreg(m_e452)


m_e51_all <-  felm(I(100*logodds(mean_xF0025))~year|country|0|year,
                   data=all_e5[
                                 (all_e5$country %in% "France" & all_e5$year==1993)==F
                               & (all_e5$country %in% "Denmark" & all_e5$year==1994)==F
                               ,])
screenreg(m_e51)

m_e451_all <-  felm(I(100*logodds(mean_xF0025))~year|country|0|year,
                    data=all_e45[
                                   (all_e45$country %in% "France" & all_e45$year==1993)==F
                                 & (all_e45$country %in% "Denmark" & all_e45$year==1994)==F
                                 ,])
screenreg(m_e451)


m_e11_all <-  felm(I(100*logodds(mean_xF0025))~year|country|0|year,
                   data=all_e1[
                                 (all_e1$country %in% "France" & all_e1$year==1993)==F
                               & (all_e1$country %in% "Denmark" & all_e1$year==1994)==F
                               ,])
screenreg(m_e11_all)


m_d0_all <-  felm(I(100*logodds(mean_xd0))~year|country|0|year,
                  data=d0[(d0$country %in% "France" & d0$year==1993)==F
                          & (d0$country %in% "Denmark" & d0$year==1994)==F
                          ,])
screenreg(m_d0)


m_beshlv_all <-  felm(I(100*logodds(pvar_between))~year|country|0|year,
                      data=all_evarlwage[(all_evarlwage$country %in% "France" & all_evarlwage$year==1993)==F
                                         & (all_evarlwage$country %in% "Denmark" & all_evarlwage$year==1994)==F
                                         ,])
screenreg(m_beshlv)


m_sh1_all <-  felm(I(100*logodds(top1share))~year|country|0|year,
                   data=all_eshare[(all_eshare$country %in% "France" & all_eshare$year==1993)==F
                                   & (all_eshare$country %in% "Denmark" & all_eshare$year==1994)==F
                                   ,])
screenreg(m_sh1)

m_sh10_all <-  felm(I(100*logodds(top10share))~year|country|0|year,
                    data=all_eshare[(all_eshare$country %in% "France" & all_eshare$year==1993)==F
                                    & (all_eshare$country %in% "Denmark" & all_eshare$year==1994)==F
                                    ,])
screenreg(m_sh10)

#TABLE 2B
screenreg(list(m_e55,m_e4545,m_e52,m_e452,m_e51,m_e451,m_e11,m_d0,m_beshlv,m_sh1,m_sh10),
          custom.model.names = c("m_e55","m_e4545","m_e52","m_e452","m_e51","m_e451","m_e11","m_d0","m_beshlv","m_sh1","m_sh10"),
          digits=1,stars=c(0.1,0.05,0.01))
htmlreg(list(m_e55,m_e4545,m_e52,m_e452,m_e51,m_e451,m_e11,m_d0,m_beshlv,m_sh1,m_sh10),
        custom.model.names = c("m_e55","m_e4545","m_e52","m_e452","m_e51","m_e451","m_e11","m_d0","m_beshlv","m_sh1","m_sh10"),
        digits=1,stars=c(0.1,0.05,0.01),file="Tables/Table 2B.html")

#TABLE 2C
screenreg(list(m_e55_all,m_e4545_all,m_e52_all,m_e452_all,m_e51_all,m_e451_all,m_e11_all,m_d0_all,m_beshlv_all,m_sh1_all,m_sh10_all),
          custom.model.names = c("m_e55","m_e4545","m_e52","m_e452","m_e51","m_e451","m_e11","m_d0","m_beshlv","m_sh1","m_sh10"),
          digits=1,stars=c(0.1,0.05,0.01))
htmlreg(list(m_e55_all,m_e4545_all,m_e52_all,m_e452_all,m_e51_all,m_e451_all,m_e11_all,m_d0_all,m_beshlv_all,m_sh1_all,m_sh10_all),
        custom.model.names = c("m_e55","m_e4545","m_e52","m_e452","m_e51","m_e451","m_e11","m_d0","m_beshlv","m_sh1","m_sh10"),
        digits=1,stars=c(0.1,0.05,0.01),file="Tables/Table 2C.html")


#------------------------------#
## Tables 3 & 4 and S8.1 & S8.2 ####
#------------------------------#


#------------------------------------------------------------------#
### Region. Global/hinterland. Data ####
#------------------------------------------------------------------#

ggc <- readRDS("Data/rds/exp_ggc.rds")

ggc$sum_w <- ggc$sum_w_xF0025
ggc$nb_obs <- ggc$N_xF0025

ggc$mean_xF0075 <- ggc$mean_xF0025+ggc$mean_xF2575
ggc$mean_xF9010 <- ggc$mean_xF9099+ggc$mean_xF9910

ggc12 <- (ggc[ggc$fwage %in% 1,c("sum_w")]*ggc[ggc$fwage %in% 1,c("mean_xF0025","mean_xF2575","mean_xF7590","mean_xF9099","mean_xF9910","mean_xF0075","mean_xF9010")]+
            ggc[ggc$fwage %in% 2,c("sum_w")]*ggc[ggc$fwage %in% 2,c("mean_xF0025","mean_xF2575","mean_xF7590","mean_xF9099","mean_xF9910","mean_xF0075","mean_xF9010")])/(ggc[ggc$fwage %in% 1,c("sum_w")]+ggc[ggc$fwage %in% 2,c("sum_w")])
ggc12$country <- ggc[ggc$fwage %in% 1,c("country")]
ggc12$year <- ggc[ggc$fwage %in% 1,c("year")]
ggc12$financial_center <- ggc[ggc$fwage %in% 1,c("financial_center")]
ggc12$sum_w <- ggc[ggc$fwage %in% 1,c("sum_w")]+ggc[ggc$fwage %in% 2,c("sum_w")]
ggc12$nb_obs <- ggc[ggc$fwage %in% 1,c("nb_obs")]+ggc[ggc$fwage %in% 2,c("nb_obs")]

ggc12$fwage <- 12

ggc45 <- (ggc[ggc$fwage %in% 4,c("sum_w")]*ggc[ggc$fwage %in% 4,c("mean_xF0025","mean_xF2575","mean_xF7590","mean_xF9099","mean_xF9910","mean_xF0075","mean_xF9010")]+
            ggc[ggc$fwage %in% 5,c("sum_w")]*ggc[ggc$fwage %in% 5,c("mean_xF0025","mean_xF2575","mean_xF7590","mean_xF9099","mean_xF9910","mean_xF0075","mean_xF9010")])/(ggc[ggc$fwage %in% 4,c("sum_w")]+ggc[ggc$fwage %in% 5,c("sum_w")])
ggc45$country <- ggc[ggc$fwage %in% 4,c("country")]
ggc45$year <- ggc[ggc$fwage %in% 4,c("year")]
ggc45$financial_center <- ggc[ggc$fwage %in% 4,c("financial_center")]
ggc45$sum_w <- ggc[ggc$fwage %in% 4,c("sum_w")]+ggc[ggc$fwage %in% 5,c("sum_w")]
ggc45$nb_obs <- ggc[ggc$fwage %in% 4,c("nb_obs")]+ggc[ggc$fwage %in% 5,c("nb_obs")]
ggc45$fwage <- 45

ggc <- rbind.fill(ggc,ggc12,ggc45)
ggc <- ggc[order(ggc$country,ggc$year,ggc$fwage,ggc$financial_center),]
rownames(ggc) <- NULL



ggc$sum_wy <- ggc$sum_w/
  ave(ggc$sum_w*(ggc$fwage %in% c(1:5))*(ggc$financial_center %in% c("",NA)==F),
      paste(ggc$country,ggc$year),FUN=function(x) sum(x,na.rm=T))

ggc$f9010_wy <- ggc$sum_w/
  ave(ggc$sum_w*(ggc$fwage %in% c(4:5))*(ggc$financial_center %in% c("",NA)==F),
      paste(ggc$country,ggc$year),
      FUN=function(x) sum(x,na.rm=T))

ggc$f9910_wy <- ggc$sum_w/
  ave(ggc$sum_w*(ggc$fwage %in% c(5))*(ggc$financial_center %in% c("",NA)==F),
      paste(ggc$country,ggc$year),FUN=function(x) sum(x,na.rm=T))


ggc$s_unit <- ave((ggc$sum_w)*(ggc$fwage %in% c(1:5)),paste(ggc$financial_center,ggc$country,ggc$year),FUN=sum)

ggc$sf9010_unit <- ave((ggc$sum_w)*(ggc$fwage %in% c(4:5)),paste(ggc$financial_center,ggc$country,ggc$year),FUN=sum)

ggc$f9010_unit <- (ggc$sf9010_unit-ifelse(ggc$fwage %in% c(4:5,45),ggc$sum_w/ggc$nb_obs,0))/(ggc$s_unit-ggc$sum_w/ggc$nb_obs)  


ggc$w_mean_xF9010_unit <- ggc$mean_xF9010-ggc$f9010_unit

ggc$w_mean_xF9010_unit_bis <- ggc$mean_xF9010-ave((ggc$sum_w)*ggc$mean_xF9010*(ggc$fwage %in% c(1:3)),paste(ggc$financial_center,ggc$country,ggc$year),FUN=sum)/
  ave((ggc$sum_w)*(ggc$fwage %in% c(1:3)),paste(ggc$financial_center,ggc$country,ggc$year),FUN=sum)

ggc$mean_xF9010_all <- ave(ggc$sum_w*ggc$mean_xF9010*(ggc$financial_center %in% c(NA,"")==F),paste(ggc$fwage,ggc$country,ggc$year),FUN=sum)/
  ave(ggc$sum_w*(ggc$financial_center %in% c(NA,"")==F),paste(ggc$fwage,ggc$country,ggc$year),FUN=sum)


ggc$f9010_unit_all <- ave(ggc$sum_w*ggc$f9010_unit*(ggc$financial_center %in% c(NA,"")==F),paste(ggc$fwage,ggc$country,ggc$year),FUN=sum)/
  ave(ggc$sum_w*(ggc$financial_center %in% c(NA,"")==F),paste(ggc$fwage,ggc$country,ggc$year),FUN=sum)

ggc$lor_f9010_unit_all <- log(oddsratio(ggc$f9010_unit_all,
                                        ave((ggc$sum_w)*(ggc$financial_center %in% c(NA,"")==T)*ggc$f9010_unit_all*(ggc$fwage %in% c(1:3)),paste(ggc$country,ggc$year),FUN=function(x) sum(x,na.rm=T))/
                                          ave((ggc$sum_w)*(ggc$financial_center %in% c(NA,"")==T)*(ggc$fwage %in% c(1:3)),paste(ggc$country,ggc$year),FUN=function(x) sum(x,na.rm=T))))

ggc$lor_w_mean_xF9010_unit <- log(oddsratio(ggc$mean_xF9010,ave((ggc$sum_w)*ggc$mean_xF9010*(ggc$fwage %in% c(1:3)),paste(ggc$financial_center,ggc$country,ggc$year),FUN=sum)/
                                              ave((ggc$sum_w)*(ggc$fwage %in% c(1:3)),paste(ggc$financial_center,ggc$country,ggc$year),FUN=sum)))

ggc$sf9910_unit <- ave((ggc$sum_w)*(ggc$fwage %in% c(5)),paste(ggc$financial_center,ggc$country,ggc$year),FUN=sum)
ggc$f9910_unit <- (ggc$sf9910_unit-ifelse(ggc$fwage %in% c(5),ggc$sum_w/ggc$nb_obs,0))/(ggc$s_unit-ggc$sum_w/ggc$nb_obs)  

ggc$w_mean_xF9910_unit <- ggc$mean_xF9910-ggc$f9910_unit

ggc$w_mean_xF9910_unit_bis <- ggc$mean_xF9910-ave((ggc$sum_w)*ggc$mean_xF9910*(ggc$fwage %in% c(1:4)),paste(ggc$financial_center,ggc$country,ggc$year),FUN=sum)/
  ave((ggc$sum_w)*(ggc$fwage %in% c(1:4)),paste(ggc$financial_center,ggc$country,ggc$year),FUN=sum)

ggc$mean_xF9910_all <- ave(ggc$sum_w*ggc$mean_xF9910*(ggc$financial_center %in% c(NA,"")==F),paste(ggc$fwage,ggc$country,ggc$year),FUN=sum)/
  ave(ggc$sum_w*(ggc$financial_center %in% c(NA,"")==F),paste(ggc$fwage,ggc$country,ggc$year),FUN=sum)


ggc$f9910_unit_all <- ave(ggc$sum_w*ggc$f9910_unit*(ggc$financial_center %in% c(NA,"")==F),paste(ggc$fwage,ggc$country,ggc$year),FUN=sum)/
  ave(ggc$sum_w*(ggc$financial_center %in% c(NA,"")==F),paste(ggc$fwage,ggc$country,ggc$year),FUN=sum)

ggc$lor_f9910_unit_all <- log(oddsratio(ggc$f9910_unit_all,
                                        ave((ggc$sum_w)*(ggc$financial_center %in% c(NA,"")==T)*ggc$f9910_unit_all*(ggc$fwage %in% c(1:4)),paste(ggc$country,ggc$year),FUN=function(x) sum(x,na.rm=T))/
                                          ave((ggc$sum_w)*(ggc$financial_center %in% c(NA,"")==T)*(ggc$fwage %in% c(1:4)),paste(ggc$country,ggc$year),FUN=function(x) sum(x,na.rm=T))))

ggc$lor_w_mean_xF9910_unit <- log(oddsratio(ggc$mean_xF9910,ave((ggc$sum_w)*ggc$mean_xF9910*(ggc$fwage %in% c(1:4)),paste(ggc$financial_center,ggc$country,ggc$year),FUN=sum)/
                                              ave((ggc$sum_w)*(ggc$fwage %in% c(1:4)),paste(ggc$financial_center,ggc$country,ggc$year),FUN=sum)))

#---------------------------------#
### Region. Global/hinterland. Table 4A & S8.2A ####
#---------------------------------#
field_all <- (ggc$financial_center  %in% c("",NA,"all"))==F & ggc$fwage %in% 45 & ggc$f9010_w>0  & is.na(ggc$f9010_w)==F &   
  (ggc$year %in% 1993 +ggc$country %in% "France")<2 &
  (ggc$year %in% 1994 +ggc$country %in% "Denmark")<2 


iie10_ggc_det_all <- felm(I(100*mean_xF9010)~factor(financial_center):year|paste(financial_center,country)|0|country+year,
                          data=ggc[field_all,],weights=ggc$f9010_w[field_all])

iiu10_ggc_det_all <- felm(I(100*f9010_unit)~factor(financial_center):year|paste(financial_center,country)|0|country+year,
                          data=ggc[field_all,],weights=ggc$f9010_w[field_all])

iie10_wggc_det_all <- felm(I(100*w_mean_xF9010_unit)~factor(financial_center):year|paste(financial_center,country)|0|country+year,
                           data=ggc[field_all,],weights=ggc$f9010_w[field_all])


screenreg(list(iie10_ggc_det_all,iiu10_ggc_det_all,iie10_wggc_det_all),
          custom.model.names=c("iie10_ggc_det","iiu10_ggc_det","iie10_wggc_det"),
          stars=c(0.1,0.05,0.01),
          digits=3)

htmlreg(list(iie10_ggc_det_all,iiu10_ggc_det_all,iie10_wggc_det_all),
        custom.model.names=c("iie10_ggc_det","iiu10_ggc_det","iie10_wggc_det"),
        stars=c(0.1,0.05,0.01),
        digits=3,
        file="Tables/Table 4A.html")

f9010_unitggc2006 <- data.frame(c("1","0")
                                ,
                                c(
                                  wtd.mean(ggc$f9010_unit[field_all & ggc$year %in% 2006 & ggc$financial_center %in% 1],weights=ggc$f9010_w[field_all  & ggc$year %in% 2006  & ggc$financial_center %in% 1]),
                                  wtd.mean(ggc$f9010_unit[field_all & ggc$year %in% 2006 & ggc$financial_center %in% 0],weights=ggc$f9010_w[field_all  & ggc$year %in% 2006  & ggc$financial_center %in% 0])))

colnames(f9010_unitggc2006) <- c("Financial_center","level_2006")
f9010_unitggc2006

f9010_unitggc2012 <- data.frame(c("1","0")
                                ,
                                c(
                                  wtd.mean(ggc$f9010_unit[field_all & ggc$year %in% 2012 & ggc$financial_center %in% 1],weights=ggc$f9010_w[field_all  & ggc$year %in% 2012  & ggc$financial_center %in% 1]),
                                  wtd.mean(ggc$f9010_unit[field_all & ggc$year %in% 2012 & ggc$financial_center %in% 0],weights=ggc$f9010_w[field_all  & ggc$year %in% 2012  & ggc$financial_center %in% 0])))

colnames(f9010_unitggc2012) <- c("Financial_center","level_2012")
f9010_unitggc2012


unlink("Tables/Table 4A-col2.html")
HTML(f9010_unitggc2012,file="Tables/Table 4A-col2.html")


#F9910
field_all <- (ggc$financial_center  %in% c("",NA,"all"))==F & ggc$fwage %in% 5 & ggc$f9910_w>0  & is.na(ggc$f9910_w)==F &   
  (ggc$year %in% 1993 +ggc$country %in% "France")<2 &
  (ggc$year %in% 1994 +ggc$country %in% "Denmark")<2 &
  (ggc$country %in% "Hungary")<1 



iie1_ggc_det_all <- felm(I(100*mean_xF9910)~factor(financial_center):year|paste(financial_center,country)|0|country+year,
                         data=ggc[field_all,],weights=ggc$f9910_w[field_all])

iiu1_ggc_det_all <- felm(I(100*f9910_unit)~factor(financial_center):year|paste(financial_center,country)|0|country+year,
                         data=ggc[field_all,],weights=ggc$f9910_w[field_all])

iie1_wggc_det_all <- felm(I(100*w_mean_xF9910_unit)~factor(financial_center):year|paste(financial_center,country)|0|country+year,
                          data=ggc[field_all,],weights=ggc$f9910_w[field_all])

# View(ggc)


screenreg(list(iie1_ggc_det_all,iiu1_ggc_det_all,iie1_wggc_det_all),
          custom.model.names=c("iie1_ggc_det","iiu1_ggc_det","iie1_wggc_det"),
          stars=c(0.1,0.05,0.01),
          digits=3)

htmlreg(list(iie1_ggc_det_all,iiu1_ggc_det_all,iie1_wggc_det_all),
        custom.model.names=c("iie1_ggc_det","iiu1_ggc_det","iie1_wggc_det"),
        stars=c(0.1,0.05,0.01),
        digits=3,
        file="Tables/Table S8.2A.html")

f9910_unitggc2006 <- data.frame(c("1","0")
                                ,
                                c(
                                  wtd.mean(ggc$f9910_unit[field_all & ggc$year %in% 2006 & ggc$financial_center %in% 1],weights=ggc$f9910_w[field_all  & ggc$year %in% 2006  & ggc$financial_center %in% 1]),
                                  wtd.mean(ggc$f9910_unit[field_all & ggc$year %in% 2006 & ggc$financial_center %in% 0],weights=ggc$f9910_w[field_all  & ggc$year %in% 2006  & ggc$financial_center %in% 0])))

colnames(f9910_unitggc2006) <- c("Financial_center","level_2006")
f9910_unitggc2006


f9910_unitggc2011 <- data.frame(c("1","0")
                                ,
                                c(
                                  wtd.mean(ggc$f9910_unit[field_all & ggc$year %in% 2011 & ggc$financial_center %in% 1],weights=ggc$f9910_w[field_all  & ggc$year %in% 2011  & ggc$financial_center %in% 1]),
                                  wtd.mean(ggc$f9910_unit[field_all & ggc$year %in% 2011 & ggc$financial_center %in% 0],weights=ggc$f9910_w[field_all  & ggc$year %in% 2011  & ggc$financial_center %in% 0])))

colnames(f9910_unitggc2011) <- c("Financial_center","level_2011")
f9910_unitggc2011

unlink("Tables/Table Table S8.2A-col2.html")
HTML(f9910_unitggc2011,file="Tables/Table S8.2A-col2.html")

#---------------------------------#
### Region. Nuts1 data ####
#---------------------------------#
nuts1 <- readRDS("Data/rds/exp_nuts1.rds")

nuts1$sum_wy <- nuts1$sum_w/
  ave(nuts1$sum_w*(nuts1$fwage %in% c(1:5))*(nuts1$nuts1 %in% c("",NA)==F),
      paste(nuts1$country,nuts1$year),FUN=function(x) sum(x,na.rm=T))

nuts1$f9010_w <- nuts1$sum_w/
  ave(nuts1$sum_w*(nuts1$fwage %in% c(4:5))*(nuts1$nuts1 %in% c("",NA)==F),
      paste(nuts1$country,nuts1$year),FUN=function(x) sum(x,na.rm=T))

nuts1$f9910_w <- nuts1$sum_w/
  ave(nuts1$sum_w*(nuts1$fwage %in% c(5))*(nuts1$nuts1 %in% c("",NA)==F),
      paste(nuts1$country,nuts1$year),FUN=function(x) sum(x,na.rm=T))

nuts1$s_unit <- ave((nuts1$sum_w)*(nuts1$fwage %in% c(1:5)),paste(nuts1$nuts1,nuts1$country,nuts1$year),FUN=sum)
nuts1$sf9010_unit <- ave((nuts1$sum_w)*(nuts1$fwage %in% c(4:5)),paste(nuts1$nuts1,nuts1$country,nuts1$year),FUN=sum)
nuts1$f9010_unit <- (nuts1$sf9010_unit-ifelse(nuts1$fwage %in% c(4:5,45),nuts1$sum_w/nuts1$nb_obs,0))/(nuts1$s_unit-nuts1$sum_w/nuts1$nb_obs)  

nuts1$w_mean_xF9010_unit <- nuts1$mean_xF9010-nuts1$f9010_unit

nuts1$w_mean_xF9010_unit_bis <- nuts1$mean_xF9010-ave((nuts1$sum_w)*nuts1$mean_xF9010*(nuts1$fwage %in% c(1:3)),paste(nuts1$nuts1,nuts1$country,nuts1$year),FUN=sum)/
  ave((nuts1$sum_w)*(nuts1$fwage %in% c(1:3)),paste(nuts1$nuts1,nuts1$country,nuts1$year),FUN=sum)

nuts1$mean_xF9010_all <- ave(nuts1$sum_w*nuts1$mean_xF9010*(nuts1$nuts1 %in% c(NA,"")==F),paste(nuts1$fwage,nuts1$country,nuts1$year),FUN=sum)/
  ave(nuts1$sum_w*(nuts1$nuts1 %in% c(NA,"")==F),paste(nuts1$fwage,nuts1$country,nuts1$year),FUN=sum)


nuts1$f9010_unit_all <- ave(nuts1$sum_w*nuts1$f9010_unit*(nuts1$nuts1 %in% c(NA,"")==F),paste(nuts1$fwage,nuts1$country,nuts1$year),FUN=sum)/
  ave(nuts1$sum_w*(nuts1$nuts1 %in% c(NA,"")==F),paste(nuts1$fwage,nuts1$country,nuts1$year),FUN=sum)

nuts1$lor_f9010_unit_all <- log(oddsratio(nuts1$f9010_unit_all,
                                          ave((nuts1$sum_w)*(nuts1$nuts1 %in% c(NA,"")==T)*nuts1$f9010_unit_all*(nuts1$fwage %in% c(1:3)),paste(nuts1$country,nuts1$year),FUN=function(x) sum(x,na.rm=T))/
                                            ave((nuts1$sum_w)*(nuts1$nuts1 %in% c(NA,"")==T)*(nuts1$fwage %in% c(1:3)),paste(nuts1$country,nuts1$year),FUN=function(x) sum(x,na.rm=T))))

nuts1$lor_w_mean_xF9010_unit <- log(oddsratio(nuts1$mean_xF9010,ave((nuts1$sum_w)*nuts1$mean_xF9010*(nuts1$fwage %in% c(1:3)),paste(nuts1$nuts1,nuts1$country,nuts1$year),FUN=sum)/
                                                ave((nuts1$sum_w)*(nuts1$fwage %in% c(1:3)),paste(nuts1$nuts1,nuts1$country,nuts1$year),FUN=sum)))



nuts1$sf9910_unit <- ave((nuts1$sum_w)*(nuts1$fwage %in% c(5)),paste(nuts1$nuts1,nuts1$country,nuts1$year),FUN=sum)
nuts1$f9910_unit <- (nuts1$sf9910_unit-ifelse(nuts1$fwage %in% 5, nuts1$sum_w/nuts1$nb_obs,0))/(nuts1$s_unit-nuts1$sum_w/nuts1$nb_obs)  

nuts1$w_mean_xF9910_unit <- nuts1$mean_xF9910-nuts1$f9910_unit

nuts1$w_mean_xF9910_unit_bis <- nuts1$mean_xF9910-ave((nuts1$sum_w)*nuts1$mean_xF9910*(nuts1$fwage %in% c(1:4)),paste(nuts1$nuts1,nuts1$country,nuts1$year),FUN=sum)/
  
  ave((nuts1$sum_w)*(nuts1$fwage %in% c(1:4)),paste(nuts1$nuts1,nuts1$country,nuts1$year),FUN=sum)

nuts1$mean_xF9910_all <- ave(nuts1$sum_w*nuts1$mean_xF9910*(nuts1$nuts1 %in% c(NA,"")==F),paste(nuts1$fwage,nuts1$country,nuts1$year),FUN=sum)/
  ave(nuts1$sum_w*(nuts1$nuts1 %in% c(NA,"")==F),paste(nuts1$fwage,nuts1$country,nuts1$year),FUN=sum)


nuts1$f9910_unit_all <- ave(nuts1$sum_w*nuts1$f9910_unit*(nuts1$nuts1 %in% c(NA,"")==F),paste(nuts1$fwage,nuts1$country,nuts1$year),FUN=sum)/
  ave(nuts1$sum_w*(nuts1$nuts1 %in% c(NA,"")==F),paste(nuts1$fwage,nuts1$country,nuts1$year),FUN=sum)

nuts1$lor_f9910_unit_all <- log(oddsratio(nuts1$f9910_unit_all,
                                          ave((nuts1$sum_w)*(nuts1$nuts1 %in% c(NA,"")==T)*nuts1$f9910_unit_all*(nuts1$fwage %in% c(1:4)),paste(nuts1$country,nuts1$year),FUN=function(x) sum(x,na.rm=T))/
                                            ave((nuts1$sum_w)*(nuts1$nuts1 %in% c(NA,"")==T)*(nuts1$fwage %in% c(1:4)),paste(nuts1$country,nuts1$year),FUN=function(x) sum(x,na.rm=T))))

nuts1$lor_w_mean_xF9910_unit <- log(oddsratio(nuts1$mean_xF9910,ave((nuts1$sum_w)*nuts1$mean_xF9910*(nuts1$fwage %in% c(1:4)),paste(nuts1$nuts1,nuts1$country,nuts1$year),FUN=sum)/
                                                ave((nuts1$sum_w)*(nuts1$fwage %in% c(1:4)),paste(nuts1$nuts1,nuts1$country,nuts1$year),FUN=sum)))


#-----------------------------#
### Region. Nuts1. Tables S8.1-2 ####
#-----------------------------#


nuts1$country <- recode(nuts1$country,
                        "Canada"="a.Canada",
                        "Denmark"="b.Denmark",
                        "Norway"="c.Norway",
                        "Sweden"="d.Sweden",
                        "France"="e.France",
                        "Netherlands"="f.Netherlands",
                        "Germany"="g.Germany",
                        "Spain"="h.Spain",
                        "Czechia"="i.Czechia",
                        "Hungary"="j.Hungary",
                        "Korea"="k.Korea",
                        "Japan"="l.Japan")

field_all <- (nuts1$nuts1  %in% c(NA,"","all",0))==F & nuts1$fwage %in% 45 & nuts1$f9010_w>0  & is.na(nuts1$f9010_w)==F &   
  (nuts1$year %in% 1993 +nuts1$country %in% "e.France")<2 &
  (nuts1$year %in% 1994 +nuts1$country %in% "b.Denmark")<2 


iie10_bb <- felm(I(100*mean_xF9010)~country:year|country|0|country+year,
                 data=nuts1[field_all,],weights=nuts1$f9010_w[field_all])


iie10_nuts1 <- felm(I(100*mean_xF9010)~country:year|paste(nuts1,country)|0|country+year,
                    data=nuts1[field_all,],weights=nuts1$f9010_w[field_all])


screenreg(list(iie10_bb,iie10_nuts1),
          custom.model.names=c("iie10_bb","iie10_nuts1"),
          stars=c(0.1,0.05,0.01),
          digits=3)

htmlreg(list(iie10_bb,iie10_nuts1),
        custom.model.names=c("iie10_bb","iie10_nuts1"),
        stars=c(0.1,0.05,0.01),
        digits=3,
        file="Tables/Table S8.1-col1-2.html")


# We drop Canada for table 3, because we don't have detailed Nuts3 estimates

field_allb <- field_all & nuts1$country %in% "a.Canada"==F

iie10_bb_allb <- felm(I(100*mean_xF9010)~year|country|0|country+year,
                      data=nuts1[field_allb,],weights=nuts1$f9010_w[field_allb])


iie10_nuts1_allb <- felm(I(100*mean_xF9010)~year|paste(nuts1,country)|0|country+year,
                         data=nuts1[field_allb,],weights=nuts1$f9010_w[field_allb])

screenreg(list(iie10_bb_allb,iie10_nuts1_allb),
          custom.model.names=c("iie10_bb","iie10_nuts1"),
          stars=c(0.1,0.05,0.01),
          digits=3)

htmlreg(list(iie10_bb_allb,iie10_nuts1_allb),
        custom.model.names=c("iie10_bb","iie10_nuts1"),
        stars=c(0.1,0.05,0.01),
        digits=3,
        file="Tables/Table 3AB-col1.html")


### Region. Nuts3. Table S8.1.3 & 3####
nuts3 <- readRDS("Data/rds/exp_nuts3.rds")

unlink("Tables/Table S8.1.col3.html")
HTML(nuts3,file="Tables/Table S8.1.col3.html")


n <- 100000
v1 <- rnorm(n,
            mean=as.numeric(gsub("\\*","",nuts3[1,2])),          
            sd=as.numeric(gsub("\\)","",gsub("\\(","",nuts3[2,2]))))
v2 <- rnorm(n,
            mean=as.numeric(gsub("*","",nuts3[3,2])),          
            sd=as.numeric(gsub("\\)","",gsub("\\(","",nuts3[4,2]))))
v3 <- rnorm(n,
            mean=as.numeric(gsub("\\*","",nuts3[5,2])),          
            sd=as.numeric(gsub("\\)","",gsub("\\(","",nuts3[6,2]))))
v4 <- rnorm(n,
            mean=as.numeric(gsub("\\*","",nuts3[7,2])),          
            sd=as.numeric(gsub("\\)","",gsub("\\(","",nuts3[8,2]))))
v5 <- rnorm(n,
            mean=as.numeric(gsub("\\*","",nuts3[9,2])),          
            sd=as.numeric(gsub("\\)","",gsub("\\(","",nuts3[10,2]))))
v6 <- rnorm(n,
            mean=as.numeric(gsub("\\*","",nuts3[11,2])),          
            sd=as.numeric(gsub("\\)","",gsub("\\(","",nuts3[12,2]))))
v7 <- rnorm(n,
            mean=as.numeric(gsub("\\*","",nuts3[13,2])),          
            sd=as.numeric(gsub("\\)","",gsub("\\(","",nuts3[14,2]))))
v8 <- rnorm(n,
            mean=as.numeric(gsub("\\*","",nuts3[15,2])),          
            sd=as.numeric(gsub("\\)","",gsub("\\(","",nuts3[16,2]))))

tapply(nuts1$f9010_w[field_all],nuts1$country[field_all],sum,na.rm=T)
colSums(table(nuts1$year[field_all],nuts1$country[field_all])>0)

v_all <- 24*v1+
  24*v2+
  29*v3+
  26*v4+
  13*v5+
  13*v6+
  15*v7+
  24*v8
v_all <- v_all/(24+24+29+26+13+13+15+24)

options(scipen=999)
mean(v_all)
sd(v_all)  
t3_C <- c(round(mean(v_all),3),paste("(",round(sd(v_all),3),")",sep=""))
unlink("Tables/Table 3C-col1.html")
HTML(data.frame(t3_C),file="Tables/Table 3C-col1.html")



### Industry 10 sectors. Data ####
sec <- readRDS("Data/rds/exp_sec.rds")


sec$sum_wy <- sec$sum_w/
  ave(sec$sum_w*(sec$fwage %in% c(1:5))*(sec$sector_agg %in% c("",NA)==F),
      paste(sec$country,sec$year),FUN=function(x) sum(x,na.rm=T))

sec$f9010_w <- sec$sum_w/
  ave(sec$sum_w*(sec$fwage %in% c(4:5))*(sec$sector_agg %in% c("",NA)==F),
      paste(sec$country,sec$year),FUN=function(x) sum(x,na.rm=T))

sec$f9910_w <- sec$sum_w/
  ave(sec$sum_w*(sec$fwage %in% c(5))*(sec$sector_agg %in% c("",NA)==F),
      paste(sec$country,sec$year),FUN=function(x) sum(x,na.rm=T))

sec$s_unit <- ave((sec$sum_w)*(sec$fwage %in% c(1:5)),paste(sec$sector_agg,sec$country,sec$year),FUN=sum)
sec$sf9010_unit <- ave((sec$sum_w)*(sec$fwage %in% c(4:5)),paste(sec$sector_agg,sec$country,sec$year),FUN=sum)
sec$f9010_unit <- (sec$sf9010_unit-ifelse(sec$fwage %in% c(4:5,45),sec$sum_w/sec$nb_obs,0))/(sec$s_unit-sec$sum_w/sec$nb_obs)  

sec$w_mean_xF9010_unit <- sec$mean_xF9010-sec$f9010_unit

sec$w_mean_xF9010_unit_bis <- sec$mean_xF9010-ave((sec$sum_w)*sec$mean_xF9010*(sec$fwage %in% c(1:3)),paste(sec$sector_agg,sec$country,sec$year),FUN=sum)/
  ave((sec$sum_w)*(sec$fwage %in% c(1:3)),paste(sec$sector_agg,sec$country,sec$year),FUN=sum)

sec$mean_xF9010_all <- ave(sec$sum_w*sec$mean_xF9010*(sec$sector_agg %in% c(NA,"")==F),paste(sec$fwage,sec$country,sec$year),FUN=sum)/
  ave(sec$sum_w*(sec$sector_agg %in% c(NA,"")==F),paste(sec$fwage,sec$country,sec$year),FUN=sum)


sec$f9010_unit_all <- ave(sec$sum_w*sec$f9010_unit*(sec$sector_agg %in% c(NA,"")==F),paste(sec$fwage,sec$country,sec$year),FUN=sum)/
  ave(sec$sum_w*(sec$sector_agg %in% c(NA,"")==F),paste(sec$fwage,sec$country,sec$year),FUN=sum)

sec$lor_f9010_unit_all <- log(oddsratio(sec$f9010_unit_all,
                                        ave((sec$sum_w)*(sec$sector_agg %in% c(NA,"")==T)*sec$f9010_unit_all*(sec$fwage %in% c(1:3)),paste(sec$country,sec$year),FUN=function(x) sum(x,na.rm=T))/
                                          ave((sec$sum_w)*(sec$sector_agg %in% c(NA,"")==T)*(sec$fwage %in% c(1:3)),paste(sec$country,sec$year),FUN=function(x) sum(x,na.rm=T))))

sec$lor_w_mean_xF9010_unit <- log(oddsratio(sec$mean_xF9010,ave((sec$sum_w)*sec$mean_xF9010*(sec$fwage %in% c(1:3)),paste(sec$sector_agg,sec$country,sec$year),FUN=sum)/
                                              ave((sec$sum_w)*(sec$fwage %in% c(1:3)),paste(sec$sector_agg,sec$country,sec$year),FUN=sum)))



sec$sf9910_unit <- ave((sec$sum_w)*(sec$fwage %in% c(5)),paste(sec$sector_agg,sec$country,sec$year),FUN=sum)
sec$f9910_unit <- (sec$sf9910_unit-ifelse(sec$fwage %in% 5, sec$sum_w/sec$nb_obs,0))/(sec$s_unit-sec$sum_w/sec$nb_obs)  

sec$w_mean_xF9910_unit <- sec$mean_xF9910-sec$f9910_unit

sec$w_mean_xF9910_unit_bis <- sec$mean_xF9910-ave((sec$sum_w)*sec$mean_xF9910*(sec$fwage %in% c(1:4)),paste(sec$sector_agg,sec$country,sec$year),FUN=sum)/
  
  ave((sec$sum_w)*(sec$fwage %in% c(1:4)),paste(sec$sector_agg,sec$country,sec$year),FUN=sum)

sec$mean_xF9910_all <- ave(sec$sum_w*sec$mean_xF9910*(sec$sector_agg %in% c(NA,"")==F),paste(sec$fwage,sec$country,sec$year),FUN=sum)/
  ave(sec$sum_w*(sec$sector_agg %in% c(NA,"")==F),paste(sec$fwage,sec$country,sec$year),FUN=sum)


sec$f9910_unit_all <- ave(sec$sum_w*sec$f9910_unit*(sec$sector_agg %in% c(NA,"")==F),paste(sec$fwage,sec$country,sec$year),FUN=sum)/
  ave(sec$sum_w*(sec$sector_agg %in% c(NA,"")==F),paste(sec$fwage,sec$country,sec$year),FUN=sum)

sec$lor_f9910_unit_all <- log(oddsratio(sec$f9910_unit_all,
                                        ave((sec$sum_w)*(sec$sector_agg %in% c(NA,"")==T)*sec$f9910_unit_all*(sec$fwage %in% c(1:4)),paste(sec$country,sec$year),FUN=function(x) sum(x,na.rm=T))/
                                          ave((sec$sum_w)*(sec$sector_agg %in% c(NA,"")==T)*(sec$fwage %in% c(1:4)),paste(sec$country,sec$year),FUN=function(x) sum(x,na.rm=T))))

sec$lor_w_mean_xF9910_unit <- log(oddsratio(sec$mean_xF9910,ave((sec$sum_w)*sec$mean_xF9910*(sec$fwage %in% c(1:4)),paste(sec$sector_agg,sec$country,sec$year),FUN=sum)/
                                              ave((sec$sum_w)*(sec$fwage %in% c(1:4)),paste(sec$sector_agg,sec$country,sec$year),FUN=sum)))

# saveRDS(sec,"Rdata/sec.rds")

# View(sec)

#---------------------------------#
### Industry 10 sectors. Tables 3, 4 & S8.1 ####
#---------------------------------#

sec$country <- recode(sec$country,
                        "Canada"="a.Canada",
                        "Denmark"="b.Denmark",
                        "Norway"="c.Norway",
                        "Sweden"="d.Sweden",
                        "France"="e.France",
                        "Netherlands"="f.Netherlands",
                        "Germany"="g.Germany",
                        "Spain"="h.Spain",
                        "Czechia"="i.Czechia",
                        "Hungary"="j.Hungary",
                        "Korea"="k.Korea",
                        "Japan"="l.Japan")
sec$sector_agg [sec$country %in% "j.Hungary" & sec$sector_ag %in% "LM" & sec$year %in% c(2003:2007)]<- "LM1"
sec$sector_agg [sec$country %in% "j.Hungary" & sec$sector_ag %in% "STU" & sec$year %in% c(2003:2007)]<- "STU1"


field_all <- ((sec$sector_agg  %in% c(NA,"","all",0)==F) &
              sec$fwage %in% 45 & sec$f9010_w>0  &
              is.na(sec$f9010_w)==F &   
             (sec$year %in% 1993 +sec$country %in% "e.France")<2 &
             (sec$year %in% 1994 +sec$country %in% "b.Denmark")<2 &
             (sec$year %in% 2017 +sec$country %in% "j.Hungary")<2 &
             (sec$country %in% c("g.Germany","i.Czechia")) <1 )
    
sec$sector_aggb <- ifelse(sec$sector_agg %in% c("K","BE","G&A","H","N"),sec$sector_agg,"ZZ")



#Trends with aggregate sectoral FE

iie10_bb_all <- felm(I(100*mean_xF9010)~year|country|0|country+year,
                     data=sec[field_all,],weights=sec$f9010_w[field_all])


iie10_sec_all <- felm(I(100*mean_xF9010)~year|paste(sector_agg,country)|0|country+year,
                      data=sec[field_all,],weights=sec$f9010_w[field_all])


screenreg(list(iie10_bb_all,iie10_sec_all),
          custom.model.names=c("iie10_bb","iie10_sec"),
          stars=c(0.1,0.05,0.01),
          digits=3)

htmlreg(list(iie10_bb_all,iie10_sec_all),
        custom.model.names=c("iie10_bb","iie10_sec"),
        stars=c(0.1,0.05,0.01),
        digits=3,
        file="Tables/Table 3AB-col2.html")

# Aggregation of detailed estimates
options(scipen=999)
X4D.sec <- readRDS("Data/rds/exp_X4D.sec.rds")

n <- 100000
v1 <- rnorm(n,
            mean=as.numeric(gsub("*","",X4D.sec[1,2])),          
            sd=as.numeric(gsub("\\)","",gsub("\\(","",X4D.sec[2,2]))))
v2 <- rnorm(n,
            mean=as.numeric(gsub("*","",X4D.sec[3,2])),          
            sd=as.numeric(gsub("\\)","",gsub("\\(","",X4D.sec[4,2]))))
v3 <- rnorm(n,
            mean=as.numeric(gsub("\\*","",X4D.sec[5,2])),          
            sd=as.numeric(gsub("\\)","",gsub("\\(","",X4D.sec[6,2]))))
v4 <- rnorm(n,
            mean=as.numeric(gsub("\\*","",X4D.sec[7,2])),          
            sd=as.numeric(gsub("\\)","",gsub("\\(","",X4D.sec[8,2]))))
v5 <- rnorm(n,
            mean=as.numeric(gsub("\\*","",X4D.sec[9,2])),          
            sd=as.numeric(gsub("\\)","",gsub("\\(","",X4D.sec[10,2]))))
v6 <- rnorm(n,
            mean=as.numeric(gsub("\\*","",X4D.sec[11,2])),          
            sd=as.numeric(gsub("\\)","",gsub("\\(","",X4D.sec[12,2]))))
v7 <- rnorm(n,
            mean=as.numeric(gsub("\\*","",X4D.sec[13,2])),          
            sd=as.numeric(gsub("\\)","",gsub("\\(","",X4D.sec[14,2]))))
v8 <- rnorm(n,
            mean=as.numeric(gsub("\\*","",X4D.sec[15,2])),          
            sd=as.numeric(gsub("\\)","",gsub("\\(","",X4D.sec[16,2]))))
v9 <- rnorm(n,
            mean=as.numeric(gsub("\\*","",X4D.sec[17,2])),          
            sd=as.numeric(gsub("\\)","",gsub("\\(","",X4D.sec[18,2]))))
v10 <- rnorm(n,
             mean=as.numeric(gsub("\\*","",X4D.sec[19,2])),          
             sd=as.numeric(gsub("\\)","",gsub("\\(","",X4D.sec[20,2]))))

tapply(sec$f9010_w[field_all],sec$country[field_all],sum,na.rm=T)
colSums(table(sec$year[field_all],sec$country[field_all])>0)

v_all <- 30*v1+
  24*v2+
  24*v3+
  29*v4+
  26*v5+
  13*v6+
  13*v7+
  14*v8+
  21*v9+
  24*v10
v_all <- v_all/(30+24+24+29+26+13+13+14+21+24)

t3_C <- c(round(mean(v_all),3),paste("(",round(sd(v_all),3),")",sep=""))

unlink("Tables/Table 3C-col2.html")
HTML(data.frame(t3_C),file="Tables/Table 3C-col2.html")


#Trends with aggregate sectoral FE per country

iie10_sec <- felm(I(100*mean_xF9010)~country:year|paste(sector_agg,country)|0|country+year,
                  data=sec[field_all,],weights=sec$f9010_w[field_all])


screenreg(list(iie10_sec),
          custom.model.names=c("iie10_sec"),
          stars=c(0.1,0.05,0.01),
          digits=3)

htmlreg(list(iie10_sec),
        custom.model.names=c("iie10_sec"),
        stars=c(0.1,0.05,0.01),
        digits=3,file="Tables/Table S8.1-col4.html")

#Trends with 4d sectoral FE per country
unlink("Tables/Table S8.1-5.html")
HTML(X4D.sec,file="Tables/Table S8.1-col5.html")


#Sectoral trends 

iie10_sec_det_all <- felm(I(100*mean_xF9010)~sector_aggb:year|paste(sector_agg,country)|0|country+year,
                          data=sec[field_all,],weights=sec$f9010_w[field_all])

iiu10_sec_det_all <- felm(I(100*f9010_unit)~sector_aggb:year|paste(sector_agg,country)|0|country+year,
                          data=sec[field_all,],weights=sec$f9010_w[field_all])

iie10_wsec_det_all <- felm(I(100*w_mean_xF9010_unit)~sector_aggb:year|paste(sector_agg,country)|0|country+year,
                           data=sec[field_all,],weights=sec$f9010_w[field_all])


screenreg(list(iie10_sec_det_all,iiu10_sec_det_all,iie10_wsec_det_all),
          custom.model.names=c("iie10_sec","iiu10_sec","iie10_wsec"),
          stars=c(0.1,0.05,0.01),
          digits=3)

htmlreg(list(iie10_sec_det_all,iiu10_sec_det_all,iie10_wsec_det_all),
        custom.model.names=c("iie10_sec","iiu10_sec","iie10_wsec"),
        stars=c(0.1,0.05,0.01),
        digits=3,
        file="Tables/Table 4B.html")


#Detail on Categorical top10% exposure 
iiu10_sec_det_allb <- felm(I(100*f9010_unit)~sector_aggb|paste(country)|0|country+year,
                           data=sec[field_all,],weights=sec$f9010_w[field_all])
summary(iiu10_sec_det_allb)

f9010_unit2006 <- data.frame(c("BE","G&A","H","K","N","ZZ")
                             ,
                             c(
                               wtd.mean(sec$f9010_unit[field_all & sec$year %in% 2006 & sec$sector_aggb %in% "BE"],weights=sec$f9010_w[field_all  & sec$year %in% 2006  & sec$sector_aggb %in% "BE"]),
                               wtd.mean(sec$f9010_unit[field_all & sec$year %in% 2006 & sec$sector_aggb %in% "G&A"],weights=sec$f9010_w[field_all  & sec$year %in% 2006  & sec$sector_aggb %in% "G&A"]),
                               wtd.mean(sec$f9010_unit[field_all & sec$year %in% 2006 & sec$sector_aggb %in% "H"],weights=sec$f9010_w[field_all  & sec$year %in% 2006  & sec$sector_aggb %in% "H"]),
                               wtd.mean(sec$f9010_unit[field_all & sec$year %in% 2006 & sec$sector_aggb %in% "K"],weights=sec$f9010_w[field_all  & sec$year %in% 2006  & sec$sector_aggb %in% "K"]),
                               wtd.mean(sec$f9010_unit[field_all & sec$year %in% 2006 & sec$sector_aggb %in% "N"],weights=sec$f9010_w[field_all  & sec$year %in% 2006  & sec$sector_aggb %in% "N"]),
                               wtd.mean(sec$f9010_unit[field_all & sec$year %in% 2006 & sec$sector_aggb %in% "ZZ"],weights=sec$f9010_w[field_all  & sec$year %in% 2006  & sec$sector_aggb %in% "ZZ"])))

colnames(f9010_unit2006) <- c("Sector","level_2006")
f9010_unit2006


f9010_unit2012 <- data.frame(c("BE","G&A","H","K","N","ZZ")
                             ,
                             c(
                               wtd.mean(sec$f9010_unit[field_all & sec$year %in% 2012 & sec$sector_aggb %in% "BE"],weights=sec$f9010_w[field_all  & sec$year %in% 2012  & sec$sector_aggb %in% "BE"]),
                               wtd.mean(sec$f9010_unit[field_all & sec$year %in% 2012 & sec$sector_aggb %in% "G&A"],weights=sec$f9010_w[field_all  & sec$year %in% 2012  & sec$sector_aggb %in% "G&A"]),
                               wtd.mean(sec$f9010_unit[field_all & sec$year %in% 2012 & sec$sector_aggb %in% "H"],weights=sec$f9010_w[field_all  & sec$year %in% 2012  & sec$sector_aggb %in% "H"]),
                               wtd.mean(sec$f9010_unit[field_all & sec$year %in% 2012 & sec$sector_aggb %in% "K"],weights=sec$f9010_w[field_all  & sec$year %in% 2012  & sec$sector_aggb %in% "K"]),
                               wtd.mean(sec$f9010_unit[field_all & sec$year %in% 2012 & sec$sector_aggb %in% "N"],weights=sec$f9010_w[field_all  & sec$year %in% 2012  & sec$sector_aggb %in% "N"]),
                               wtd.mean(sec$f9010_unit[field_all & sec$year %in% 2012 & sec$sector_aggb %in% "ZZ"],weights=sec$f9010_w[field_all  & sec$year %in% 2012  & sec$sector_aggb %in% "ZZ"])))

colnames(f9010_unit2012) <- c("Sector","level_2012")
f9010_unit2012

unlink("Tables/Table 4B-col2.html")
HTML(f9010_unit2012,file="Tables/Table 4B-col2.html")



#Top 1%
field_all <- (sec$sector_agg  %in% c(NA,"","all",0))==F & 
  sec$fwage %in% 5 & sec$f9910_w>0  & 
  is.na(sec$f9910_w)==F &   
  (sec$year %in% 1993 +sec$country %in% "e.France")<2 &
  (sec$year %in% 1994 +sec$country %in% "b.Denmark")<2 &
  (sec$year %in% 2017 +sec$country %in% "j.Hungary")<2 &
  (sec$country %in% c("g.Germany","i.Czechia")+0)<1 


iie1_sec_det_all <- felm(I(100*mean_xF9910)~sector_aggb:year|paste(sector_agg,country)|0|country+year,
                          data=sec[field_all,],weights=sec$f9910_w[field_all])

iiu1_sec_det_all <- felm(I(100*f9910_unit)~sector_aggb:year|paste(sector_agg,country)|0|country+year,
                          data=sec[field_all,],weights=sec$f9910_w[field_all])

iie1_wsec_det_all <- felm(I(100*w_mean_xF9910_unit)~sector_aggb:year|paste(sector_agg,country)|0|country+year,
                           data=sec[field_all,],weights=sec$f9910_w[field_all])


screenreg(list(iie1_sec_det_all,iiu1_sec_det_all,iie1_wsec_det_all),
          custom.model.names=c("iie1_sec_det","iiu1_sec_det","iie1_wsec_det"),
          stars=c(0.1,0.05,0.01),
          digits=3)

htmlreg(list(iie1_sec_det_all,iiu1_sec_det_all,iie1_wsec_det_all),
        custom.model.names=c("iie1_sec","iiu1_sec","iie1_wsec"),
        stars=c(0.1,0.05,0.01),
        digits=3,
        file="Tables/Table S8.2B.html")



#Detail on Categorical top1% exposure 
f9910_unit2006 <- data.frame(c("BE","G&A","H","K","N","ZZ")
                             ,
                             c(
wtd.mean(sec$f9910_unit[field_all & sec$year %in% 2006 & sec$sector_aggb %in% "BE"],weights=sec$f9910_w[field_all  & sec$year %in% 2006  & sec$sector_aggb %in% "BE"]),
wtd.mean(sec$f9910_unit[field_all & sec$year %in% 2006 & sec$sector_aggb %in% "G&A"],weights=sec$f9910_w[field_all  & sec$year %in% 2006  & sec$sector_aggb %in% "G&A"]),
wtd.mean(sec$f9910_unit[field_all & sec$year %in% 2006 & sec$sector_aggb %in% "H"],weights=sec$f9910_w[field_all  & sec$year %in% 2006  & sec$sector_aggb %in% "H"]),
wtd.mean(sec$f9910_unit[field_all & sec$year %in% 2006 & sec$sector_aggb %in% "K"],weights=sec$f9910_w[field_all  & sec$year %in% 2006  & sec$sector_aggb %in% "K"]),
wtd.mean(sec$f9910_unit[field_all & sec$year %in% 2006 & sec$sector_aggb %in% "N"],weights=sec$f9910_w[field_all  & sec$year %in% 2006  & sec$sector_aggb %in% "N"]),
wtd.mean(sec$f9910_unit[field_all & sec$year %in% 2006 & sec$sector_aggb %in% "ZZ"],weights=sec$f9910_w[field_all  & sec$year %in% 2006  & sec$sector_aggb %in% "ZZ"])))

colnames(f9910_unit2006) <- c("Sector","level_2006")
f9910_unit2006


f9910_unit2011 <- data.frame(c("BE","G&A","H","K","N","ZZ")
                             ,
                             c(
                               wtd.mean(sec$f9910_unit[field_all & sec$year %in% 2012 & sec$sector_aggb %in% "BE"],weights=sec$f9910_w[field_all  & sec$year %in% 2011  & sec$sector_aggb %in% "BE"]),
                               wtd.mean(sec$f9910_unit[field_all & sec$year %in% 2011 & sec$sector_aggb %in% "G&A"],weights=sec$f9910_w[field_all  & sec$year %in% 2011  & sec$sector_aggb %in% "G&A"]),
                               wtd.mean(sec$f9910_unit[field_all & sec$year %in% 2011 & sec$sector_aggb %in% "H"],weights=sec$f9910_w[field_all  & sec$year %in% 2011  & sec$sector_aggb %in% "H"]),
                               wtd.mean(sec$f9910_unit[field_all & sec$year %in% 2011 & sec$sector_aggb %in% "K"],weights=sec$f9910_w[field_all  & sec$year %in% 2011  & sec$sector_aggb %in% "K"]),
                               wtd.mean(sec$f9910_unit[field_all & sec$year %in% 2011 & sec$sector_aggb %in% "N"],weights=sec$f9910_w[field_all  & sec$year %in% 2011  & sec$sector_aggb %in% "N"]),
                               wtd.mean(sec$f9910_unit[field_all & sec$year %in% 2011 & sec$sector_aggb %in% "ZZ"],weights=sec$f9910_w[field_all  & sec$year %in% 2011  & sec$sector_aggb %in% "ZZ"])))

colnames(f9910_unit2011) <- c("Sector","level_2011")
f9910_unit2011

unlink("Tables/Table S8.2B-col2.html")
HTML(f9910_unit2011,file="Tables/Table S8.2B-col2.html")



#-------------------------------#
## Table 5. Size decrease #####
#-------------------------------#


estfe <- readRDS("Data/rds/exp_estfe.rds")



estfe <- data.frame(sapply(estfe,FUN=function(x) gsub("\\)","",gsub("\\(","",x))))
estfe <- data.frame(sapply(estfe,FUN=function(x) gsub("\\*","",x)))
estfe0 <-estfe

estfe0[,c(2:15)] <- sapply(estfe0[,c(2:15)],FUN=function(x) 100*as.numeric(x))

estfe[,c(2:15)] <- sapply(estfe[,c(2:15)],FUN=function(x) round(100*as.numeric(x),3))
estfe[,c(3,5,7,9,11,13,15)] <- sapply(estfe[,c(3,5,7,9,11,13,15)],FUN=function(x) paste("(",x,")",sep=""))

estfe
estfe_star <- estfe0[,c(2,4,6,8,10,12,14)]/estfe0[,c(3,5,7,9,11,13,15)]
estfe_star2 <- sapply(estfe_star,function(x) ifelse(abs(x)>2.576,"***",ifelse(abs(x)>1.96,"**",ifelse(abs(x)>1.645,"*",""))))

estfe[,2] <- paste(estfe[,2],estfe_star2[,1], sep="")
estfe[,4] <- paste(estfe[,4],estfe_star2[,2], sep="")
estfe[,6] <- paste(estfe[,6],estfe_star2[,3], sep="")
estfe[,8] <- paste(estfe[,8],estfe_star2[,4], sep="")
estfe[,10] <- paste(estfe[,10],estfe_star2[,5], sep="")
estfe[,12] <- paste(estfe[,12],estfe_star2[,6], sep="")
estfe[,14] <- paste(estfe[,14],estfe_star2[,7], sep="")

testfe <-t(estfe)

colnames(testfe) <- c(
  "Canada",
  "Denmark",
  "Norway",
  "Sweden",
  "France",
  "Netherlands",
  "Spain",
  "Hungary",
  "Korea",
  "Japan")
testfe

estfe2 <- estfe[,c("Country","estfesz_sz","estfesz_sz_se",
                   "estfesz2_sz","estfesz2_sz_se","estfesz2_sz_neg","estfesz2_sz_neg_se")]

estfe2$parm <- paste(estfe2$estfesz_sz,estfe2$estfesz_sz_se)
estfe2$parm2 <- paste(estfe2$estfesz2_sz,estfe2$estfesz2_sz_se)
estfe2$parm3 <- paste(estfe2$estfesz2_sz_neg,estfe2$estfesz2_sz_neg_se)
  estfe3 <- estfe2[,c("Country","parm","parm2","parm3")]

unlink("Tables/Table 5.html")
HTML(estfe3,file="Tables/Table 5.html")



# TABLE 5B
n <- 1000000
v1 <- rnorm(n,
            mean=estfe0[1,8],          
            sd=estfe0[1,9])
v2 <- rnorm(n,
            mean=estfe0[2,8],          
            sd=estfe0[2,9])
v3 <- rnorm(n,
            mean=estfe0[3,8],          
            sd=estfe0[3,9])
v4 <- rnorm(n,
            mean=estfe0[4,8],          
            sd=estfe0[4,9])
v5 <- rnorm(n,
            mean=estfe0[5,8],          
            sd=estfe0[5,9])
v6 <- rnorm(n,
            mean=estfe0[6,8],          
            sd=estfe0[6,9])
v7 <- rnorm(n,
            mean=estfe0[7,8],          
            sd=estfe0[7,9])
v8 <- rnorm(n,
            mean=estfe0[8,8],          
            sd=estfe0[8,9])
v9 <- rnorm(n,
            mean=estfe0[9,8],          
            sd=estfe0[9,9])
v10 <- rnorm(n,
             mean=estfe0[10,8],          
             sd=estfe0[10,9])


v_all <- 30*v1+
  24*v2+
  24*v3+
  29*v4+
  26*v5+
  13*v6+
  13*v7+
  16*v8+
  21*v9+
  24*v10
v_all <- v_all/(30+24+24+29+26+13+13+16
                +21+24)

mean(v_all)
sd(v_all)  
t_5B1 <- c(round(mean(v_all),3),paste("(",round(sd(v_all),3),")",sep=""))
t_5B1

n <- 1000000
v1 <- rnorm(n,
            mean=estfe0[1,12],          
            sd=estfe0[1,13])
v2 <- rnorm(n,
            mean=estfe0[2,12],          
            sd=estfe0[2,13])
v3 <- rnorm(n,
            mean=estfe0[3,12],          
            sd=estfe0[3,13])
v4 <- rnorm(n,
            mean=estfe0[4,12],          
            sd=estfe0[4,13])
v5 <- rnorm(n,
            mean=estfe0[5,12],          
            sd=estfe0[5,13])
v6 <- rnorm(n,
            mean=estfe0[6,12],          
            sd=estfe0[6,13])
v7 <- rnorm(n,
            mean=estfe0[7,12],          
            sd=estfe0[7,13])
v8 <- rnorm(n,
            mean=estfe0[8,12],          
            sd=estfe0[8,13])
v9 <- rnorm(n,
            mean=estfe0[9,12],          
            sd=estfe0[9,13])
v10 <- rnorm(n,
             mean=estfe0[10,12],          
             sd=estfe0[10,13])


v_all <- 30*v1+
  24*v2+
  24*v3+
  29*v4+
  26*v5+
  13*v6+
  13*v7+
  16*v8+
  21*v9+
  24*v10
v_all <- v_all/(30+24+24+29+26+13+13+16
                +21+24)

t_5B2 <- c(round(mean(v_all),3),paste("(",round(sd(v_all),3),")",sep=""))
t_5B2


n <- 1000000
v1 <- rnorm(n,
            mean=estfe0[1,14],          
            sd=estfe0[1,15])
v2 <- rnorm(n,
            mean=estfe0[2,14],          
            sd=estfe0[2,15])
v3 <- rnorm(n,
            mean=estfe0[3,14],          
            sd=estfe0[3,15])
v4 <- rnorm(n,
            mean=estfe0[4,14],          
            sd=estfe0[4,15])
v5 <- rnorm(n,
            mean=estfe0[5,14],          
            sd=estfe0[5,15])
v6 <- rnorm(n,
            mean=estfe0[6,14],          
            sd=estfe0[6,15])
v7 <- rnorm(n,
            mean=estfe0[7,14],          
            sd=estfe0[7,15])
v8 <- rnorm(n,
            mean=estfe0[8,14],          
            sd=estfe0[8,15])
v9 <- rnorm(n,
            mean=estfe0[9,14],          
            sd=estfe0[9,15])
v10 <- rnorm(n,
             mean=estfe0[10,14],          
             sd=estfe0[10,15])


v_all <- 30*v1+
  24*v2+
  24*v3+
  29*v4+
  26*v5+
  13*v6+
  13*v7+
  16*v8+
  21*v9+
  24*v10
v_all <- v_all/(30+24+24+29+26+13+13+16
                +21+24)

t_5B3 <- c(round(mean(v_all),3),paste("(",round(sd(v_all),3),")",sep=""))
t_5B3

table_5B <- data.frame(t_5B1,t_5B2,t_5B3)
unlink("Tables/Table 5B.html")
HTML(table_5B,file="Tables/Table 5B.html")


## Table 3. Last column  & Table S8-6 ####

#Benchmark
my_f <- 
  (all_e45$country %in% "France" & all_e45$year==1993)==F & 
  (all_e45$country %in% "Denmark" & all_e45$year==1994)==F&
  all_e45$country %in% c("Germany","Czechia")==F

m_e4545<- felm(I(100*(mean_xF9010))~year|country|0|year,
               data=all_e45[my_f,])
screenreg(m_e4545,digits=3)
htmlreg(m_e4545,digits=3,file="Tables/Table 3A-col3.html")


# With workplace fixed effect. Aggregate effect
estfe0
n <- 1000000
v1 <- rnorm(n,
            mean=estfe0[1,4],          
            sd=estfe0[1,5])
v2 <- rnorm(n,
            mean=estfe0[2,4],          
            sd=estfe0[2,5])
v3 <- rnorm(n,
            mean=estfe0[3,4],          
            sd=estfe0[3,5])
v4 <- rnorm(n,
            mean=estfe0[4,4],          
            sd=estfe0[4,5])
v5 <- rnorm(n,
            mean=estfe0[5,4],          
            sd=estfe0[5,5])
v6 <- rnorm(n,
            mean=estfe0[6,4],          
            sd=estfe0[6,5])
v7 <- rnorm(n,
            mean=estfe0[7,4],          
            sd=estfe0[7,5])
v8 <- rnorm(n,
            mean=estfe0[8,4],          
            sd=estfe0[8,5])
v9 <- rnorm(n,
            mean=estfe0[9,4],          
            sd=estfe0[9,5])
v10 <- rnorm(n,
             mean=estfe0[10,4],          
             sd=estfe0[10,5])


v_all <- 30*v1+
  24*v2+
  24*v3+
  29*v4+
  26*v5+
  13*v6+
  13*v7+
  16*v8+
  21*v9+
  24*v10
v_all <- v_all/(30+24+24+29+26+13+13+16
                +21+24)

t_3C3 <- c(round(mean(v_all),3),paste("(",round(sd(v_all),3),")",sep=""))
t_3C3

options(scipen=999)
unlink("Tables/Table 3C-col3.html")
HTML(t_3C3,file="Tables/Table 3C-col3.html")

# TABLE S8.1-6
estfe3 <- rbind.fill(estfe[,c("Country","estfe")],data.frame(estfe[,c("estfe_se")]))
estfe3$id <- ifelse(as.numeric(row.names(estfe3))<11,as.numeric(row.names(estfe3)),as.numeric(row.names(estfe3))-10)

estfe3$parm <- ifelse(as.numeric(row.names(estfe3))<11,estfe3[,2],estfe3[,3])
estfe3 <- estfe3[order(estfe3$id),c("Country","parm")]

unlink("Tables/Table S8.1-col6.html")
HTML(estfe3,file="Tables/Table S8.1-col6.html")

#---------------------------------#
#---------------------------------#
## Table 6. Overall Regression ####
#---------------------------------#
#---------------------------------#

pp <- readRDS("Data/rds/exp_pp.rds")

pp$lf9010xf9010 <- simplelag(pp$f9010xf9010,pp$country)
pp$dvolgdp <- pp$volgdp - simplelag(pp$volgdp,pp$country)
pp$dvolgdp[is.na(pp$dvolgdp)==T] <-0
pp$dvolgdp_cumneg <- ave(pp$dvolgdp*(pp$dvolgdp<0),pp$country,FUN=function(x) cumsum(x))


# Dependent variable
pp$f9010xf9010_dm <- pp$f9010xf9010-ave(pp$f9010xf9010,pp$country,FUN=function(x) mean(x,na.rm=T))
pp$f9010xf9010_std <- pp$f9010xf9010_dm/sd(pp$f9010xf9010_dm,na.rm=T)
pp$f9010xf9010_100 <- 100*pp$f9010xf9010


# Deindustrialisation indicator
pp$ma_empshare_dm <- pp$ma_empshare-ave(pp$ma_empshare,pp$country,FUN=function(x) mean(x,na.rm=T))
pp$ma_empshare_std <- pp$ma_empshare_dm/sd(pp$ma_empshare_dm,na.rm=T)

# Finance indicator
pp$volgdp_dm <- pp$volgdp-ave(pp$volgdp,pp$country,FUN=function(x) mean(x,na.rm=T))
pp$volgdp_std <- pp$volgdp_dm/sd(pp$volgdp_dm,na.rm=T)


# Global city indicator
pp$gc_wshare_dm <- pp$gc_wshare-ave(pp$gc_wshare,pp$country,FUN=function(x) mean(x,na.rm=T))
pp$gc_wshare_std <- pp$gc_wshare_dm/sd(pp$gc_wshare_dm,na.rm=T)

pp$gc_wshare3 <- ifelse(pp$country %in% "Korea",pp$Wage_Sh_Seoul,pp$gc_wshare)
pp$gc_wshare3_dm <- pp$gc_wshare3-ave(pp$gc_wshare3,pp$country,FUN=function(x) mean(x,na.rm=T))
pp$gc_wshare3_std <- pp$gc_wshare3_dm/sd(pp$gc_wshare3_dm,na.rm=T)


#Average size indicator
pp$lm_orgsize_dm <- log(pp$m_orgsize)-ave(log(pp$m_orgsize),pp$country,FUN=function(x) mean(x,na.rm=T))
pp$lm_orgsize_std <- pp$lm_orgsize_dm/sd(pp$lm_orgsize_dm,na.rm=T)

pp$lm_orgsize_std_neg <- ifelse(pp$lm_orgsize_std-simplelag(pp$lm_orgsize_std,pp$country)<0,
                                pp$lm_orgsize_std-simplelag(pp$lm_orgsize_std,pp$country),0)
pp$lm_orgsize_std_neg <- ifelse(is.na(pp$lm_orgsize_std_neg),0,pp$lm_orgsize_std_neg)
pp$lm_orgsize_std_cumneg <- ave(pp$lm_orgsize_std_neg,pp$country,FUN=cumsum)

# ICT
pp$ict_share_asset2 <- ifelse(is.na(pp$ict_share_asset_oecd),pp$ict_share_asset,pp$ict_share_asset_oecd)
pp$ict_share_asset_source <- 1-is.na(pp$ict_share_asset_oecd)

pp$ict_share_asset2_dm <- pp$ict_share_asset2-ave(pp$ict_share_asset2,pp$country,FUN=function(x) mean(x,na.rm=T))
pp$ict_share_asset2_std <- pp$ict_share_asset2_dm/sd(pp$ict_share_asset2_dm,na.rm=T)


# FDI
pp$fdi_dm <- pp$fdi-ave(pp$fdi,pp$country,FUN=function(x) mean(x,na.rm=T))
pp$fdi_std <- pp$fdi_dm/sd(pp$fdi_dm,na.rm=T)


#Log wage OECD
pp$loecd_wage_dm <- log(pp$oecd_wage)-ave(log(pp$oecd_wage),pp$country,FUN=function(x) mean(x,na.rm=T))
pp$loecd_wage_std <- pp$loecd_wage_dm/sd(pp$loecd_wage_dm,na.rm=T)

#Log population OECD
pp$lpop2064_dm <- log(pp$pop2064)-ave(log(pp$pop2064),pp$country,FUN=function(x) mean(x,na.rm=T))
pp$lpop2064_std <- pp$lpop2064_dm/sd(pp$lpop2064_dm,na.rm=T)

# Inequality
pp$var_total_dm <- pp$var_total-ave(pp$var_total,pp$country,FUN=function(x) mean(x,na.rm=T))
pp$var_total_std <- pp$var_total_dm/sd(pp$var_total_dm,na.rm=T)

pp$var_between_dm <- pp$var_between-ave(pp$var_between,pp$country,FUN=function(x) mean(x,na.rm=T))
pp$var_between_std <- pp$var_between_dm/sd(pp$var_between_dm,na.rm=T)

pp$var_within_dm <- pp$var_within-ave(pp$var_within,pp$country,FUN=function(x) mean(x,na.rm=T))
pp$var_within_std <- pp$var_within_dm/sd(pp$var_within_dm,na.rm=T)

pp$pvar_between_dm <- pp$pvar_between-ave(pp$pvar_between,pp$country,FUN=function(x) mean(x,na.rm=T))
pp$pvar_between_std <- pp$pvar_between_dm/sd(pp$pvar_between_dm,na.rm=T)

pp$top10share_dm <- pp$top10share-ave(pp$top10share,pp$country,FUN=function(x) mean(x,na.rm=T))
pp$top10share_std <- pp$top10share_dm/sd(pp$top10share_dm,na.rm=T)


# EXPLANATIONS
ex00 <- felm(f9010xf9010_std~
               #Deindustrialization
               # + ma_empshare_std
             
             #Financialization
             # + volgdp_std
             
             #Size
             # + lm_orgsize_std
             # + lm_orgsize_std_cumneg
             
             #Global cities
             # + gc_wshare3_std             
             
             # ICT share assets
             # + ict_share_asset2_std
             
             # FDI
             # + fdi_std
             
             # Controls
             + loecd_wage_std
             + lpop2064_std
             
             |country+year|0|year,data=pp[,])


ex01 <- felm(f9010xf9010_std~
               #Deindustrialization
               + ma_empshare_std
             
             #Financialization
             # + volgdp_std
             
             #Size
             # + lm_orgsize_std
             # + lm_orgsize_std_cumneg
             
             #Global cities
             # + gc_wshare3_std             
             
             # ICT share assets
              # + ict_share_asset2_std
             
             # FDI
             # + fdi_std
             
             # Controls
             + loecd_wage_std
             + lpop2064_std
             
             |country+year|0|year,data=pp[,])

screenreg(ex01)
ex02 <- felm(f9010xf9010_std~
               #Deindustrialization
               # + ma_empshare_std
               
               #Financialization
             + volgdp_std

             #Size
             # + lm_orgsize_std
             # + lm_orgsize_std_cumneg
              
             #Global cities
             # + gc_wshare3_std
             
             # ICT share assets
              # + ict_share_asset2_std
             
             # FDI
             # + fdi_std

             # Controls
             + loecd_wage_std
             + lpop2064_std
             |country+year|0|year,data=pp[,])
screenreg(ex02)

ex03 <- felm(f9010xf9010_std~
               #Deindustrialization
               # + ma_empshare_std
               
               #Financialization
             # + volgdp_std

             #Size
             + lm_orgsize_std

             #Global cities
             # + gc_wshare3_std
             
             # ICT share assets
              # + ict_share_asset2_std
             
             # FDI
             # + fdi_std
             
             # Controls
             + loecd_wage_std
             + lpop2064_std
             |country+year|0|year,data=pp[,])
screenreg(ex03,digits=3)

ex04 <- felm(f9010xf9010_std~
          #Deindustrialization
               # + ma_empshare_std
               
               #Financialization
             # + volgdp_std
             
             #Size
             + lm_orgsize_std
             + lm_orgsize_std_cumneg
             
             #Global cities
             # + gc_wshare3_std
             
             # ICT share assets
              # + ict_share_asset2_std
             
             # FDI
             # + fdi_std
             
             # Controls
             + loecd_wage_std
             + lpop2064_std
             |country+year|0|year,data=pp[,])
summary(ex04)

ex05 <- felm(f9010xf9010_std~
               #YEAR LINEAR TREND
               # + year 
               
               #Deindustrialization
               # + ma_empshare_std
               
               #Financialization
             # + volgdp_std
             
             #Size
              # + lm_orgsize_std
             # + lm_orgsize_std_cumneg
             
             
             #Global cities
             + gc_wshare3_std             
             
             # ICT share assets
              # + ict_share_asset2_std
             
             # FDI
             # + fdi_std
             
             # Controls
             + loecd_wage_std
             + lpop2064_std
             |country+year|0|year,data=pp[,])

ex06 <- felm(f9010xf9010_std~
               #Deindustrialization
               # + ma_empshare_std
               
               #Financialization
             # + volgdp_std
             
             #Size
               # + lm_orgsize_std
               # + lm_orgsize_std_cumneg
               
             
             #Global cities
             # + gc_wshare3_std
             
             # ICT share assets
              # + ict_share_asset2_std
             
             # FDI
             + fdi_std
             
             # Controls
             + loecd_wage_std
             + lpop2064_std
             |country+year|0|year,data=pp[,])

ex07 <- felm(f9010xf9010_std~
               #Deindustrialization
               # + ma_empshare_std
               
               #Financialization
             # + volgdp_std
             
             #Size
               # + lm_orgsize_std
               # + lm_orgsize_std_cumneg
               
             
             #Global cities
             # + gc_wshare3_std
             
             # ICT share assets
             + ict_share_asset2_std

             # FDI
             # + fdi_std
             
             # Controls
             + loecd_wage_std
             + lpop2064_std
             |country+year|0|year,data=pp[,])

screenreg(ex07)

screenreg(list(ex00,ex01,ex03,ex04,ex06,ex05,ex07,ex02),stars=c(0.1,0.05,0.01),digits=3)
htmlreg(list(ex00,ex01,ex03,ex04,ex06,ex05,ex07,ex02),stars=c(0.1,0.05,0.01),digits=3,file="Tables/Table 6A.html")

ex11 <- felm(f9010xf9010_std~
               #Deindustrialization
               + ma_empshare_std
             
             #Financialization
             # + volgdp_std
             
             #Size
             + lm_orgsize_std
             + lm_orgsize_std_cumneg
             
             #Global cities
             # + gc_wshare3_std
             
             # FDI
             + fdi_std
             
             # ICT share assets
              # + ict_share_asset2_std
             
             
             # Controls
             + loecd_wage_std
             + lpop2064_std
             
             |country+year|0|year,data=pp[,])

ex12 <- felm(f9010xf9010_std~
               #Deindustrialization
               + ma_empshare_std
             
             #Financialization
             # + volgdp_std
             
             #Size
             + lm_orgsize_std
             + lm_orgsize_std_cumneg
             
             
             #Global cities
             # + gc_wshare3_std
             
             # ICT share assets
              # + ict_share_asset2_std
             
             # FDI
             + fdi_std
             
             # Controls
             + loecd_wage_std
             + lpop2064_std
             |country+year|0|year,data=pp[is.na(pp$gc_wshare3_std)==F,])

ex13 <- felm(f9010xf9010_std~
               #Deindustrialization
               + ma_empshare_std
             
             #Financialization
             # + volgdp_std
             
             #Size
             + lm_orgsize_std
             + lm_orgsize_std_cumneg
             
             #Global cities
             + gc_wshare3_std
             
             # ICT share assets
              # + ict_share_asset2_std
             
             # FDI
             + fdi_std

             # Controls
             + loecd_wage_std
             + lpop2064_std
             |country+year|0|year,data=pp[,])

ex14 <- felm(f9010xf9010_std~
               #Deindustrialization
               + ma_empshare_std
             
             #Financialization
             # + volgdp_std
             
             #Size
             + lm_orgsize_std
             + lm_orgsize_std_cumneg
             
             
             #Global cities
             + gc_wshare3_std
             
             # ICT share assets
              # + ict_share_asset2_std
             
             # FDI
             + fdi_std
             
             # Controls
             + loecd_wage_std
             + lpop2064_std
             |country+year|0|year,data=pp[is.na(pp$ict_share_asset2_std)==F,])

ex15 <- felm(f9010xf9010_std~
               #Deindustrialization
               + ma_empshare_std

             #Financialization
             # + volgdp_std
             
             #Size
             + lm_orgsize_std
             + lm_orgsize_std_cumneg
             
             
             #Global cities
             + gc_wshare3_std
             
             # ICT share assets
             + ict_share_asset2_std
             
             # FDI
             + fdi_std
             
             # Controls
             + loecd_wage_std
             + lpop2064_std
             |country+year|0|year,data=pp[,])
screenreg(ex15,stars=c(0.1,0.05,0.01),digits=3)

ex16 <- felm(f9010xf9010_std~
               #Deindustrialization
               + ma_empshare_std

             #Financialization
             # + volgdp_std
             
             #Size
             + lm_orgsize_std
             + lm_orgsize_std_cumneg
             
             
             #Global cities
             + gc_wshare3_std
             
             # ICT share assets
             + ict_share_asset2_std
             
             # FDI
             + fdi_std
             
             # Controls
             + loecd_wage_std
             + lpop2064_std
             |country+year|0|year,data=pp[is.na(pp$volgdp_std)==F,])


ex17 <- felm(f9010xf9010_std~
               #Deindustrialization
               + ma_empshare_std

             #Financialization
             + volgdp_std
             
             #Size
             + lm_orgsize_std
             + lm_orgsize_std_cumneg
             
             #Global cities
             + gc_wshare3_std
             
             # ICT share assets
             + ict_share_asset2_std
             
             # FDI
             + fdi_std

             # Controls
             + loecd_wage_std
             + lpop2064_std
             |country+year|0|year,data=pp[,])
summary(ex17)
screenreg(list(ex00,ex11,ex12,ex13,ex14,ex15,ex16,ex17),stars=c(0.1,0.05,0.01),digits=3)
htmlreg(list(ex00,ex11,ex12,ex13,ex14,ex15,ex16,ex17),stars=c(0.1,0.05,0.01),digits=3,file="Tables/Table 6B.html")

#-------------------------##
## Table S6. Exposure and inequality  ####
#-------------------------##
in01 <- felm(f9010xf9010_std~top10share_std|country+year|0|year,data=pp[,])
in02 <- felm(f9010xf9010_std~var_total_std |country+year|0|year,data=pp[,])
in03 <- felm(f9010xf9010_std~var_within_std|country+year|0|year,data=pp[,])
in04 <- felm(f9010xf9010_std~var_between_std |country+year|0|year,data=pp[,])
in05 <- felm(f9010xf9010_std~pvar_between_std |country+year|0|year,data=pp[,])


screenreg(list(in01,in02,in03,in04,in05),stars=c(0.1,0.05,0.01),digits=3)
htmlreg(list(in01,in02,in03,in04,in05),stars=c(0.1,0.05,0.01),digits=3,file="Tables/Table S6.html")


#------------------------------------------------------------------#
## Table S4.1 Robustness ####
#------------------------------------------------------------------#
all_e <- readRDS("Data/rds/exp_all_e.rds")
all_f <- readRDS("Data/rds/exp_all_f.rds")
all_h <- readRDS("Data/rds/exp_all_h.rds")

all_e45 <- all_e[all_e$fwage==45,]
all_f45 <- all_f[all_f$fwage==45,]
all_h45 <- all_h[all_h$fhwage==45,]

all_e45$country <- recode(all_e45$country,
                      "Canada"="a.Canada",
                      "Denmark"="b.Denmark",
                      "Norway"="c.Norway",
                      "Sweden"="d.Sweden",
                      "France"="e.France",
                      "Netherlands"="f.Netherlands",
                      "Germany"="g.Germany",
                      "Spain"="h.Spain",
                      "Czechia"="i.Czechia",
                      "Hungary"="j.Hungary",
                      "Korea"="k.Korea",
                      "Japan"="l.Japan")

all_f45$country <- recode(all_f45$country,
                          "Canada"="a.Canada",
                          "Denmark"="b.Denmark",
                          "Norway"="c.Norway",
                          "Sweden"="d.Sweden",
                          "France"="e.France",
                          "Netherlands"="f.Netherlands",
                          "Germany"="g.Germany",
                          "Spain"="h.Spain",
                          "Czechia"="i.Czechia",
                          "Hungary"="j.Hungary",
                          "Korea"="k.Korea",
                          "Japan"="l.Japan")

all_h45$country <- recode(all_h45$country,
                          "Canada"="a.Canada",
                          "Denmark"="b.Denmark",
                          "Norway"="c.Norway",
                          "Sweden"="d.Sweden",
                          "France"="e.France",
                          "Netherlands"="f.Netherlands",
                          "Germany"="g.Germany",
                          "Spain"="h.Spain",
                          "Czechia"="i.Czechia",
                          "Hungary"="j.Hungary",
                          "Korea"="k.Korea",
                          "Japan"="l.Japan")

m_e4545 <-  felm(I(100*logodds(mean_xF9010))~year:country|country|0|year,
                 data=all_e45[(all_e45$country %in% "France" & all_e45$year==1993)==F
                              & (all_e45$country %in% "Denmark" & all_e45$year==1994)==F
                              ,])
screenreg(m_e4545)


m_e451 <-  felm(I(100*logodds(mean_xF0025))~year:country|country|0|year,
                data=all_e45[(all_e45$country %in% "France" & all_e45$year==1993)==F
                             & (all_e45$country %in% "Denmark" & all_e45$year==1994)==F
                             ,])
screenreg(m_e451)


m_f4545 <-  felm(I(100*logodds(mean_xF9010))~year:country|country|0|year,
                 data=all_f45[(all_f45$country %in% "France" & all_f45$year==1993)==F
                              & (all_f45$country %in% "Denmark" & all_f45$year==1994)==F
                              ,])
screenreg(m_f4545)


m_f451 <-  felm(I(100*logodds(mean_xF0025))~year:country|country|0|year,
                data=all_f45[(all_f45$country %in% "France" & all_f45$year==1993)==F
                             & (all_f45$country %in% "Denmark" & all_f45$year==1994)==F
                             ,])
screenreg(m_f451)


m_h4545 <-  felm(I(100*logodds(mean_xF9010h))~year:country|country|0|year,
                 data=all_h45[(all_h45$country %in% "France" & all_h45$year==1993)==F
                              & (all_h45$country %in% "Denmark" & all_h45$year==1994)==F
                              ,])
screenreg(m_h4545)


m_h451 <-  felm(I(100*logodds(mean_xF0025h))~year:country|country|0|year,
                data=all_h45[(all_h45$country %in% "France" & all_h45$year==1993)==F
                             & (all_h45$country %in% "Denmark" & all_h45$year==1994)==F
                             ,])
screenreg(m_h451)


m_e4545_all <-  felm(I(100*logodds(mean_xF9010))~year|country|0|year,
                     data=all_e45[(all_e45$country %in% "France" & all_e45$year==1993)==F
                                  & (all_e45$country %in% "Denmark" & all_e45$year==1994)==F
                                  ,])
screenreg(m_e4545_all)


m_e451_all <-  felm(I(100*logodds(mean_xF0025))~year|country|0|year,
                    data=all_e45[(all_e45$country %in% "France" & all_e45$year==1993)==F
                                 & (all_e45$country %in% "Denmark" & all_e45$year==1994)==F
                                 ,])
screenreg(m_e451_all)


m_f4545_all <-  felm(I(100*logodds(mean_xF9010))~year|country|0|year,
                     data=all_f45[(all_f45$country %in% "France" & all_f45$year==1993)==F
                                  & (all_f45$country %in% "Denmark" & all_f45$year==1994)==F
                                  ,])
screenreg(m_f4545_all)


m_f451_all <-  felm(I(100*logodds(mean_xF0025))~year|country|0|year,
                    data=all_f45[(all_f45$country %in% "France" & all_f45$year==1993)==F
                                 & (all_f45$country %in% "Denmark" & all_f45$year==1994)==F
                                 ,])
screenreg(m_f451_all)


m_h4545_all <-  felm(I(100*logodds(mean_xF9010h))~year|country|0|year,
                     data=all_h45[(all_h45$country %in% "France" & all_h45$year==1993)==F
                                  & (all_h45$country %in% "Denmark" & all_h45$year==1994)==F
                                  ,])
screenreg(m_h4545_all)


m_h451_all <-  felm(I(100*logodds(mean_xF0025h))~year|country|0|year,
                    data=all_h45[(all_h45$country %in% "France" & all_h45$year==1993)==F
                                 & (all_h45$country %in% "Denmark" & all_h45$year==1994)==F
                                 ,])
screenreg(m_h451_all)

#Table S4.1A
screenreg(list(m_e4545,m_f4545,m_h4545,m_e451,m_f451,m_h451),
          digits=1,stars=c(0.1,0.05,0.01))
htmlreg(list(m_e4545,m_f4545,m_h4545,m_e451,m_f451,m_h451),
        digits=1,stars=c(0.1,0.05,0.01),file="Tables/Table S4.1A.html")

#Table S4.1b
screenreg(list(m_e4545_all,m_f4545_all,m_h4545_all,m_e451_all,m_f451_all,m_h451_all),
          digits=1,stars=c(0.1,0.05,0.01))
htmlreg(list(m_e4545_all,m_f4545_all,m_h4545_all,m_e451_all,m_f451_all,m_h451_all),
        digits=1,stars=c(0.1,0.05,0.01),file="Tables/Table S4.1B.html")

#---------------------------------#
## Table S5.  Specificity ####
#---------------------------------#

### Data preparation ####

# Earnings relative exposure 
all_e <- readRDS("Data/rds/exp_all_e.rds")

eo <- all_e


colnames(eo)[which(names(eo) %in% c("mean_xF0025","mean_xF2575","mean_xF7590","mean_xF9099","mean_xF9910","mean_xF0075","mean_xF9010"))] <- c("xF0025","xF2575","xF7590","xF9099","xF9910","xF0075","xF9010")

eo_p <- eo[eo$fwage %in% NA,c("year","country","xF0025","xF2575","xF7590","xF9099","xF9910","xF0075","xF9010")]

colnames(eo_p)[which(names(eo_p) %in% c("xF0025","xF2575","xF7590","xF9099","xF9910","xF0075","xF9010"))] <- c("F0025","F2575","F7590","F9099","F9910","F0075","F9010")

eo_p$t <- eo_p$F0025+eo_p$F2575+eo_p$F7590+eo_p$F9099+eo_p$F9910
eo <- merge(eo,eo_p,by=c("year","country"),all.x=TRUE)
str(eo)

eo_w1 <- eo[eo$fwage %in% 1,c("year","country","sum_w")]
eo_w2 <- eo[eo$fwage %in% 2,c("sum_w")]
eo_w3 <- eo[eo$fwage %in% 3,c("sum_w")]
eo_w4 <- eo[eo$fwage %in% 4,c("sum_w")]
eo_w5 <- eo[eo$fwage %in% 5,c("sum_w")]

eo_w12 <- eo[eo$fwage %in% 12,c("year","sum_w")]

eo_wall <- eo[eo$fwage %in% NA,c("sum_w")]

rownames(eo_w1) <- NULL

eo_w <- cbind(eo_w1,eo_w2,eo_w3,eo_w4,eo_w5,eo_wall)  

colnames(eo_w) <- c("year","country","weotF0025","weotF2575","weotF7590","weotF9099","weotF9910","sweot")  

eo_w$weotF0025 <- eo_w$weotF0025/eo_w$sweot
eo_w$weotF2575 <- eo_w$weotF2575/eo_w$sweot
eo_w$weotF7590 <- eo_w$weotF7590/eo_w$sweot
eo_w$weotF9099 <- eo_w$weotF9099/eo_w$sweot
eo_w$weotF9910 <- eo_w$weotF9910/eo_w$sweot
eo_w$weotF0075 <- eo_w$weotF0025+eo_w$weotF2575
eo_w$weotF9010 <- eo_w$weotF9910+eo_w$weotF9099

eo <- merge(eo,eo_w,by=c("year","country"),all.x=TRUE)

eo$xF0025_n <- ifelse(eo$fwage==1,(eo$F0025-eo$weotF0025*eo$xF0025)/(1-eo$weotF0025),
                      ifelse(eo$fwage==2,(eo$F0025-eo$weotF2575*eo$xF0025)/(1-eo$weotF2575),
                             ifelse(eo$fwage==3,(eo$F0025-eo$weotF7590*eo$xF0025)/(1-eo$weotF7590),
                                    ifelse(eo$fwage==4,(eo$F0025-eo$weotF9099*eo$xF0025)/(1-eo$weotF9099),
                                           ifelse(eo$fwage==5,(eo$F0025-eo$weotF9910*eo$xF0025)/(1-eo$weotF9910),
                                                  ifelse(eo$fwage==12,(eo$F0025-eo$weotF0075*eo$xF0025)/(1-eo$weotF0075),
                                                         ifelse(eo$fwage==45,(eo$F0025-eo$weotF9010*eo$xF0025)/(1-eo$weotF9010),NA)))))))

eo$xF2575_n <- ifelse(eo$fwage==1,(eo$F2575-eo$weotF0025*eo$xF2575)/(1-eo$weotF0025),
                      ifelse(eo$fwage==2,(eo$F2575-eo$weotF2575*eo$xF2575)/(1-eo$weotF2575),
                             ifelse(eo$fwage==3,(eo$F2575-eo$weotF7590*eo$xF2575)/(1-eo$weotF7590),
                                    ifelse(eo$fwage==4,(eo$F2575-eo$weotF9099*eo$xF2575)/(1-eo$weotF9099),
                                           ifelse(eo$fwage==5,(eo$F2575-eo$weotF9910*eo$xF2575)/(1-eo$weotF9910),
                                                  ifelse(eo$fwage==12,(eo$F2575-eo$weotF0075*eo$xF2575)/(1-eo$weotF0075),
                                                         ifelse(eo$fwage==45,(eo$F2575-eo$weotF9010*eo$xF2575)/(1-eo$weotF9010),NA)))))))


eo$xF7590_n <- ifelse(eo$fwage==1,(eo$F7590-eo$weotF0025*eo$xF7590)/(1-eo$weotF0025),
                      ifelse(eo$fwage==2,(eo$F7590-eo$weotF2575*eo$xF7590)/(1-eo$weotF2575),
                             ifelse(eo$fwage==3,(eo$F7590-eo$weotF7590*eo$xF7590)/(1-eo$weotF7590),
                                    ifelse(eo$fwage==4,(eo$F7590-eo$weotF9099*eo$xF7590)/(1-eo$weotF9099),
                                           ifelse(eo$fwage==5,(eo$F7590-eo$weotF9910*eo$xF7590)/(1-eo$weotF9910),
                                                  ifelse(eo$fwage==12,(eo$F7590-eo$weotF0075*eo$xF7590)/(1-eo$weotF0075),
                                                         ifelse(eo$fwage==45,(eo$F7590-eo$weotF9010*eo$xF7590)/(1-eo$weotF9010),NA)))))))

eo$xF9099_n <- ifelse(eo$fwage==1,(eo$F9099-eo$weotF0025*eo$xF9099)/(1-eo$weotF0025),
                      ifelse(eo$fwage==2,(eo$F9099-eo$weotF2575*eo$xF9099)/(1-eo$weotF2575),
                             ifelse(eo$fwage==3,(eo$F9099-eo$weotF7590*eo$xF9099)/(1-eo$weotF7590),
                                    ifelse(eo$fwage==4,(eo$F9099-eo$weotF9099*eo$xF9099)/(1-eo$weotF9099),
                                           ifelse(eo$fwage==5,(eo$F9099-eo$weotF9910*eo$xF9099)/(1-eo$weotF9910),
                                                  ifelse(eo$fwage==12,(eo$F9099-eo$weotF0075*eo$xF9099)/(1-eo$weotF0075),
                                                         ifelse(eo$fwage==45,(eo$F9099-eo$weotF9010*eo$xF9099)/(1-eo$weotF9010),NA)))))))


eo$xF9910_n <- ifelse(eo$fwage==1,(eo$F9910-eo$weotF0025*eo$xF9910)/(1-eo$weotF0025),
                      ifelse(eo$fwage==2,(eo$F9910-eo$weotF2575*eo$xF9910)/(1-eo$weotF2575),
                             ifelse(eo$fwage==3,(eo$F9910-eo$weotF7590*eo$xF9910)/(1-eo$weotF7590),
                                    ifelse(eo$fwage==4,(eo$F9910-eo$weotF9099*eo$xF9910)/(1-eo$weotF9099),
                                           ifelse(eo$fwage==5,(eo$F9910-eo$weotF9910*eo$xF9910)/(1-eo$weotF9910),
                                                  ifelse(eo$fwage==12,(eo$F9910-eo$weotF0075*eo$xF9910)/(1-eo$weotF0075),
                                                         ifelse(eo$fwage==45,(eo$F9910-eo$weotF9010*eo$xF9910)/(1-eo$weotF9010),NA)))))))


eo$xF0075_n <- ifelse(eo$fwage==1,(eo$F0075-eo$weotF0025*eo$xF0075)/(1-eo$weotF0025),
                      ifelse(eo$fwage==2,(eo$F0075-eo$weotF2575*eo$xF0075)/(1-eo$weotF2575),
                             ifelse(eo$fwage==3,(eo$F0075-eo$weotF7590*eo$xF0075)/(1-eo$weotF7590),
                                    ifelse(eo$fwage==4,(eo$F0075-eo$weotF9099*eo$xF0075)/(1-eo$weotF9099),
                                           ifelse(eo$fwage==5,(eo$F0075-eo$weotF9910*eo$xF0075)/(1-eo$weotF9910),
                                                  ifelse(eo$fwage==12,(eo$F0075-eo$weotF0075*eo$xF0075)/(1-eo$weotF0075),
                                                         ifelse(eo$fwage==45,(eo$F0075-eo$weotF9010*eo$xF0075)/(1-eo$weotF9010),NA)))))))

eo$xF9010_n <- ifelse(eo$fwage==1,(eo$F9010-eo$weotF0025*eo$xF9010)/(1-eo$weotF0025),
                      ifelse(eo$fwage==2,(eo$F9010-eo$weotF2575*eo$xF9010)/(1-eo$weotF2575),
                             ifelse(eo$fwage==3,(eo$F9010-eo$weotF7590*eo$xF9010)/(1-eo$weotF7590),
                                    ifelse(eo$fwage==4,(eo$F9010-eo$weotF9099*eo$xF9010)/(1-eo$weotF9099),
                                           ifelse(eo$fwage==5,(eo$F9010-eo$weotF9910*eo$xF9010)/(1-eo$weotF9010),
                                                  ifelse(eo$fwage==12,(eo$F9010-eo$weotF0075*eo$xF9010)/(1-eo$weotF0075),
                                                         ifelse(eo$fwage==45,(eo$F9010-eo$weotF9010*eo$xF9010)/(1-eo$weotF9010),NA)))))))


eo$odF0025 <- ((eo$xF0025/(1-eo$xF0025))/((eo$xF0025_n/(1-eo$xF0025_n))))
eo$odF2575 <- ((eo$xF2575/(1-eo$xF2575))/((eo$xF2575_n/(1-eo$xF2575_n))))
eo$odF7590 <- ((eo$xF7590/(1-eo$xF7590))/((eo$xF7590_n/(1-eo$xF7590_n))))
eo$odF9099 <- ((eo$xF9099/(1-eo$xF9099))/((eo$xF9099_n/(1-eo$xF9099_n))))
eo$odF9910 <- ((eo$xF9910/(1-eo$xF9910))/((eo$xF9910_n/(1-eo$xF9910_n))))
eo$odF0075 <- ((eo$xF0075/(1-eo$xF0075))/((eo$xF0075_n/(1-eo$xF0075_n))))
eo$odF9010 <- ((eo$xF9010/(1-eo$xF9010))/((eo$xF9010_n/(1-eo$xF9010_n))))


eo$lodF0025 <- log(eo$odF0025)
eo$lodF2575 <- log(eo$odF2575)
eo$lodF7590 <- log(eo$odF7590)
eo$lodF9099 <- log(eo$odF9099)
eo$lodF9910 <- log(eo$odF9910)
eo$lodF0075 <- log(eo$odF0075)
eo$lodF9010 <- log(eo$odF9010)


eo$dF0025 <- (eo$xF0025-eo$xF0025_n)
eo$dF2575 <- (eo$xF2575-eo$xF2575_n)
eo$dF7590 <- (eo$xF7590-eo$xF7590_n)
eo$dF9099 <- (eo$xF9099-eo$xF9099_n)
eo$dF9910 <- (eo$xF9910-eo$xF9910_n)
eo$dF0075 <- (eo$xF0075-eo$xF0075_n)
eo$dF9010 <- (eo$xF9010-eo$xF9010_n)

all_eo1 <- eo[eo$fwage %in% 1 & eo$year %in% c(1980:1989,2001.1)==F,]
all_eo45 <- eo[eo$fwage %in% 45  & eo$year %in% c(1980:1989,2001.1)==F,]
all_eo5 <- eo[eo$fwage %in% 5  & eo$year %in% c(1980:1989,2001.1)==F,]

# Data migration
seg_mig <- readRDS("Data/rds/exp_seg_mig.rds")

seg_mig1 <- seg_mig[seg_mig$migrant %in% 1,]
seg_mig0 <- seg_mig[seg_mig$migrant %in% 0,]

seg_mig1$omigrant <- oddsratio(seg_mig1$mean_xmigrant,seg_mig0$mean_xmigrant)
seg_mig1$lomigrant <- log(seg_mig1$omigrant)

# Data gender
seg_gdr <- readRDS("Data/rds/exp_seg_gdr.rds")

seg_gdr <- seg_gdr[(seg_gdr$country=="Denmark")+(seg_gdr$year==1994)<2,]
seg_gdr1 <- seg_gdr[seg_gdr$female %in% 1,]
seg_gdr0 <- seg_gdr[seg_gdr$female %in% 0,]
seg_gdr1$ofemale <- oddsratio(seg_gdr1$mean_xfemale,seg_gdr0$mean_xfemale)
seg_gdr1$lofemale <- log(seg_gdr1$ofemale)

# Data Age
age <- readRDS("Data/rds/exp_seg_age.rds")


colnames(age)[which(names(age) %in% c("mean_xage1","mean_xage2","mean_xage3","mean_xage4"))] <- c("xage1","xage2","xage3","xage4")
age$xage_t <- age$xage1+age$xage2+age$xage3+age$xage4
age$xage1 <- age$xage1/age$xage_t
age$xage2 <- age$xage2/age$xage_t
age$xage3 <- age$xage3/age$xage_t
age$xage4 <- age$xage4/age$xage_t
age$xage_t <- age$xage1+age$xage2+age$xage3+age$xage4  

age_p <- age[age$age %in% NA,c("year","country","xage1","xage2","xage3","xage4")]
colnames(age_p)[which(names(age_p) %in% c("xage1","xage2","xage3","xage4"))] <- c("age1","age2","age3","age4")
age_p$t <- age_p$age1+age_p$age2+age_p$age3+age_p$age4
age <- merge(age,age_p,by=c("year","country"),all.x=TRUE)

age_w1 <- age[age$age %in% 1,c("year","country","sum_w_xage1")]
age_w2 <- age[age$age %in% 2,c("year","country","sum_w_xage1")]
age_w3 <- age[age$age %in% 3,c("year","country","sum_w_xage1")]
age_w4 <- age[age$age %in% 4,c("year","country","sum_w_xage1")]
age_wall <- age[age$age %in% NA,c("year","country","sum_w_xage1")]

age_w <- merge(age_w1,age_w2,by=c("year","country"),all=TRUE)
colnames(age_w) <- c("year","country","wgtage1","wgtage2")
age_w <- merge(age_w,age_w3,by=c("year","country"),all=TRUE)
colnames(age_w) <- c("year","country","wgtage1","wgtage2","wgtage3")
age_w <- merge(age_w,age_w4,by=c("year","country"),all=TRUE)
colnames(age_w) <- c("year","country","wgtage1","wgtage2","wgtage3","wgtage4")
age_w <- merge(age_w,age_wall,by=c("year","country"),all=TRUE)
colnames(age_w) <- c("year","country","wgtage1","wgtage2","wgtage3","wgtage4","swgt")  

age_w$wgtage1 <- age_w$wgtage1/age_w$swgt
age_w$wgtage2 <- age_w$wgtage2/age_w$swgt
age_w$wgtage3 <- age_w$wgtage3/age_w$swgt
age_w$wgtage4 <- age_w$wgtage4/age_w$swgt

age <- merge(age,age_w,by=c("year","country"),all.x=TRUE)

str(age)
age$xage1_n <- ifelse(age$age==1,(age$age1-age$wgtage1*age$xage1)/(1-age$wgtage1),
                      ifelse(age$age==2,(age$age1-age$wgtage2*age$xage1)/(1-age$wgtage2),
                             ifelse(age$age==3,(age$age1-age$wgtage3*age$xage1)/(1-age$wgtage3),
                                    ifelse(age$age==4,(age$age1-age$wgtage4*age$xage1)/(1-age$wgtage4),NA))))

age$xage2_n <- ifelse(age$age==1,(age$age2-age$wgtage1*age$xage2)/(1-age$wgtage1),
                      ifelse(age$age==2,(age$age2-age$wgtage2*age$xage2)/(1-age$wgtage2),
                             ifelse(age$age==3,(age$age2-age$wgtage3*age$xage2)/(1-age$wgtage3),
                                    ifelse(age$age==4,(age$age2-age$wgtage4*age$xage2)/(1-age$wgtage4),NA))))


age$xage3_n <- ifelse(age$age==1,(age$age3-age$wgtage1*age$xage3)/(1-age$wgtage1),
                      ifelse(age$age==2,(age$age3-age$wgtage2*age$xage3)/(1-age$wgtage2),
                             ifelse(age$age==3,(age$age3-age$wgtage3*age$xage3)/(1-age$wgtage3),
                                    ifelse(age$age==4,(age$age3-age$wgtage4*age$xage3)/(1-age$wgtage4),NA))))                                  

age$xage4_n <- ifelse(age$age==1,(age$age4-age$wgtage1*age$xage4)/(1-age$wgtage1),
                      ifelse(age$age==2,(age$age4-age$wgtage2*age$xage4)/(1-age$wgtage2),
                             ifelse(age$age==3,(age$age4-age$wgtage3*age$xage4)/(1-age$wgtage3),
                                    ifelse(age$age==4,(age$age4-age$wgtage4*age$xage4)/(1-age$wgtage4),NA))))                                  




age$odage1 <- ((age$xage1/(1-age$xage1))/((age$xage1_n/(1-age$xage1_n))))
age$odage2 <- ((age$xage2/(1-age$xage2))/((age$xage2_n/(1-age$xage2_n))))
age$odage3 <- ((age$xage3/(1-age$xage3))/((age$xage3_n/(1-age$xage3_n))))
age$odage4 <- ((age$xage4/(1-age$xage4))/((age$xage4_n/(1-age$xage4_n))))


age$lodage1 <- log(age$odage1)
age$lodage2 <- log(age$odage2)
age$lodage3 <- log(age$odage3)
age$lodage4 <- log(age$odage4)


age$dage1 <- (age$xage1-age$xage1_n)
age$dage2 <- (age$xage2-age$xage2_n)
age$dage3 <- (age$xage3-age$xage3_n)
age$dage4 <- (age$xage4-age$xage4_n)


seg_agena <- age[is.na(age$age),]
seg_age4 <- age[age$age==4,]
seg_age3 <- age[age$age==3,]
seg_age1 <- age[age$age==1,]

# Data occupation
occ <- readRDS("Data/rds/exp_seg_occ.rds")

colnames(occ)[which(names(occ) %in% c("mean_xocc1","mean_xocc2","mean_xocc3"))] <- c("xocc1","xocc2","xocc3")

occ_p <- occ[occ$occ %in% NA,c("year","country","xocc1","xocc2","xocc3")]
colnames(occ_p)[which(names(occ_p) %in% c("xocc1","xocc2","xocc3"))] <- c("occ1","occ2","occ3")
occ_p$t <- occ_p$occ1+occ_p$occ2+occ_p$occ3
occ <- merge(occ,occ_p,by=c("year","country"),all.x=TRUE)

occ_w1 <- occ[occ$occ %in% 1,c("year","country","sum_w_xocc1")]
occ_w2 <- occ[occ$occ %in% 2,c("sum_w_xocc1")]
occ_w3 <- occ[occ$occ %in% 3,c("sum_w_xocc1")]
occ_wall <- occ[is.na(occ$occ),c("sum_w_xocc1")]

rownames(occ_w1) <- NULL
occ_w <- cbind(occ_w1,occ_w2,occ_w3,occ_wall)  
colnames(occ_w) <- c("year","country","wgtocc1","wgtocc2","wgtocc3","swgt")  

occ_w$wgtocc1 <- occ_w$wgtocc1/occ_w$swgt
occ_w$wgtocc2 <- occ_w$wgtocc2/occ_w$swgt
occ_w$wgtocc3 <- occ_w$wgtocc3/occ_w$swgt

occ <- merge(occ,occ_w,by=c("year","country"),all.x=TRUE)

occ$xocc1_n <- ifelse(occ$occ==1,(occ$occ1-occ$wgtocc1*occ$xocc1)/(1-occ$wgtocc1),
                      ifelse(occ$occ==2,(occ$occ1-occ$wgtocc2*occ$xocc1)/(1-occ$wgtocc2),
                             ifelse(occ$occ==3,(occ$occ1-occ$wgtocc3*occ$xocc1)/(1-occ$wgtocc3),NA)))

occ$xocc2_n <- ifelse(occ$occ==1,(occ$occ2-occ$wgtocc1*occ$xocc2)/(1-occ$wgtocc1),
                      ifelse(occ$occ==2,(occ$occ2-occ$wgtocc2*occ$xocc2)/(1-occ$wgtocc2),
                             ifelse(occ$occ==3,(occ$occ2-occ$wgtocc3*occ$xocc2)/(1-occ$wgtocc3),NA)))

occ$xocc3_n <- ifelse(occ$occ==1,(occ$occ3-occ$wgtocc1*occ$xocc3)/(1-occ$wgtocc1),
                      ifelse(occ$occ==2,(occ$occ3-occ$wgtocc2*occ$xocc3)/(1-occ$wgtocc2),
                             ifelse(occ$occ==3,(occ$occ3-occ$wgtocc3*occ$xocc3)/(1-occ$wgtocc3),NA)))


occ$odocc1 <- ((occ$xocc1/(1-occ$xocc1))/((occ$xocc1_n/(1-occ$xocc1_n))))
occ$odocc2 <- ((occ$xocc2/(1-occ$xocc2))/((occ$xocc2_n/(1-occ$xocc2_n))))
occ$odocc3 <- ((occ$xocc3/(1-occ$xocc3))/((occ$xocc3_n/(1-occ$xocc3_n))))


occ$lodocc1 <- log(occ$odocc1)
occ$lodocc2 <- log(occ$odocc2)
occ$lodocc3 <- log(occ$odocc3)

occ$docc1 <- (occ$xocc1-occ$xocc1_n)
occ$docc2 <- (occ$xocc2-occ$xocc2_n)
occ$docc3 <- (occ$xocc3-occ$xocc3_n)


seg_occna <- occ[is.na(occ$occ),]

seg_occ1 <- occ[occ$occ %in% 1
                & ((occ$country=="Denmark")+(occ$year==1994 | occ$year>2009))<2
                & ((occ$country=="Germany")+(occ$year>=2011))<2
                & ((occ$country=="Korea")+(occ$year>=2008))<2
                & ((occ$country=="Korea")+(occ$year<=1992))<2
                ,]
seg_occ3 <- occ[occ$occ %in% 3
                & ((occ$country=="Denmark")+(occ$year==1994 | occ$year>2009))<2
                & ((occ$country=="Germany")+(occ$year>=2011))<2
                & ((occ$country=="Korea")+(occ$year>=2008))<2
                & ((occ$country=="Korea")+(occ$year<=1992))<2
                ,]

### Trend curves ####
g_eo4545 <- mygraph2(all_eo45$lodF9010,"Top 10%  establishment relative exposure to top 10%","Relative exposure (odds-ratio)","",ratio="logodds",legpos="bottomright",myyaxt="n")
g_eo55 <- mygraph2(all_eo5$lodF9910,"Top 1%  establishment relative exposure to top 1% ","Relative exposure (odds-ratio)","",ratio="logodds",legpos="bottomright",myyaxt="n")
mig_eo <- mygraph2(seg_mig1$lomigrant,"Migrant establishment isolation",
                   "Relative exposure (odds-ratio)","",
                   ratio="logodds",legpos="topright",
                   myyaxt="n")
gdr_eo <- mygraph2(seg_gdr1$lofemale,"Female establishment isolation",
                   "Relative exposure (odds-ratio)","",
                   ratio="logodds",legpos="topright",
                   myyaxt="n")

g_oage4 <- mygraph2(seg_age4$lodage4,"Age >55 establishment isolation","Relative establishement exposure  (odds-ratio)","",ratio="logodds",
                    legpos="topright",myyaxt="n")

g_oage1 <- mygraph2(seg_age1$lodage1,"Age <30 establishment isolation","Relative establishement exposure  (odds-ratio)","",ratio="logodds",
                    legpos="topleft",myyaxt="n")

g_oocc1 <- mygraph2(seg_occ1$lodocc1,"Managers and professionals relative establishment isolation",
                    "Relative establishement exposure  (odds-ratio)","",ratio="logodds",legpos="topleft",myyaxt="n")

g_oocc3 <- mygraph2(seg_occ3$lodocc3,"Working-class employees relative establishment isolation",
                    "Relative establishement exposure  (odds-ratio)","",ratio="logodds",legpos="topleft",myyaxt="n")


### Table S5 ####
table_S5 <- rbind.fill(data.frame(t(g_eo55$des[-1,c("last","gr_rate")])),
                      data.frame(t(g_eo4545$des[-1,c("last","gr_rate")])),
                      data.frame(t(mig_eo$des[-1,c("last","gr_rate")])),
                      data.frame(t(gdr_eo$des[-1,c("last","gr_rate")])),
                      data.frame(t(g_oage4$des[-1,c("last","gr_rate")])),
                      data.frame(t(g_oage1$des[-1,c("last","gr_rate")])),
                      # data.frame(t(g_oeduc$des[-1,c("last","gr_rate")])),
                      data.frame(t(g_oocc1$des[-1,c("last","gr_rate")])),
                      data.frame(t(g_oocc3$des[-1,c("last","gr_rate")])))

colnames(table_S5) <- g_eo55$des$country[-1]
table_S5$type <- rep(c("last","gr_rate"),8)

table_S5_l <- exp(table_S5[table_S5$type %in% "last",c(1:13)])
table_S5_e <- table_S5[table_S5$type %in% "gr_rate",c(1:13)]
table_S5_l$seg <- c("Top 1%","Top 10%","Migrant","Female","Age>55","Age<31","Mgrs & Pro.","Execution workers")
table_S5_e$seg <- c("Top 1%","Top 10%","Migrant","Female","Age>55","Age<31","Mgrs & Pro.","Execution workers")


table_S5A <- data.frame(t(table_S5_l[,c("seg","Canada","Denmark","Norway","Sweden","France","Netherlands","Germany","Spain",
                           "Czechia","Hungary","Korea","Japan","Adj. Mean")]))


#table S5A Export
write.csv(table_S5A,"Tables/Table S5A.csv",row.names = T, na=".")

#table 5 Regressions
m_eo55 <-  felm(I(100*lodF9910)~year:country|country|0|year,
                data=all_eo5[(all_eo5$country %in% "France" & all_eo5$year==1993)==F
                             & (all_eo5$country %in% "Denmark" & all_eo5$year==1994)==F
                             ,])

m_eo4545 <-  felm(I(100*lodF9010)~year:country|country|0|year,
                  data=all_eo45[(all_eo45$country %in% "France" & all_eo45$year==1993)==F
                                & (all_eo45$country %in% "Denmark" & all_eo45$year==1994)==F
                                ,])
summary(m_eo4545)

m_mig_eo <-  felm(I(100*lomigrant)~year:country|country|0|year,
                  data=seg_mig1[(seg_mig1$country %in% "France" & seg_mig1$year==1993)==F
                                & (seg_mig1$country %in% "Denmark" & seg_mig1$year==1994)==F
                                ,])

m_gdr_eo <-  felm(I(100*lofemale)~year:country|country|0|year,
                  data=seg_gdr1[(seg_gdr1$country %in% "France" & seg_gdr1$year==1993)==F
                                & (seg_gdr1$country %in% "Denmark" & seg_gdr1$year==1994)==F
                                ,])

m_oage1 <-  felm(I(100*lodage1)~year:country|country|0|year,
                 data=seg_age1[(seg_age1$country %in% "France" & seg_age1$year==1993)==F
                               & (seg_age1$country %in% "Denmark" & seg_age1$year==1994)==F
                               ,])

m_oage4 <-  felm(I(100*lodage4)~year:country|country|0|year,
                 data=seg_age4[(seg_age4$country %in% "France" & seg_age4$year==1993)==F
                               & (seg_age4$country %in% "Denmark" & seg_age4$year==1994)==F
                               ,])

m_oocc1 <-  felm(I(100*lodocc1)~year:country|country|0|year,
                 data=seg_occ1[(seg_occ1$country %in% "France" & seg_occ1$year==1993)==F
                               & (seg_occ1$country %in% "Denmark" & seg_occ1$year==1994)==F
                               ,])

m_oocc3 <-  felm(I(100*lodocc3)~year:country|country|0|year,
                 data=seg_occ3[(seg_occ3$country %in% "France" & seg_occ3$year==1993)==F
                               & (seg_occ3$country %in% "Denmark" & seg_occ3$year==1994)==F
                               ,])

m_eo55_all <-  felm(I(100*lodF9910)~year|country|0|year,
                    data=all_eo5[(all_eo5$country %in% "France" & all_eo5$year==1993)==F
                                 & (all_eo5$country %in% "Denmark" & all_eo5$year==1994)==F
                                 ,])

m_eo4545_all <-  felm(I(100*lodF9010)~year|country|0|year,
                      data=all_eo45[(all_eo45$country %in% "France" & all_eo45$year==1993)==F
                                    & (all_eo45$country %in% "Denmark" & all_eo45$year==1994)==F
                                    ,])

m_mig_eo_all <-  felm(I(100*lomigrant)~year|country|0|year,
                      data=seg_mig1[(seg_mig1$country %in% "France" & seg_mig1$year==1993)==F
                                    & (seg_mig1$country %in% "Denmark" & seg_mig1$year==1994)==F
                                    ,])

m_gdr_eo_all <-  felm(I(100*lofemale)~year|country|0|year,
                      data=seg_gdr1[(seg_gdr1$country %in% "France" & seg_gdr1$year==1993)==F
                                    & (seg_gdr1$country %in% "Denmark" & seg_gdr1$year==1994)==F
                                    ,])

m_oage1_all <-  felm(I(100*lodage1)~year|country|0|year,
                     data=seg_age1[(seg_age1$country %in% "France" & seg_age1$year==1993)==F
                                   & (seg_age1$country %in% "Denmark" & seg_age1$year==1994)==F
                                   ,])

m_oage4_all <-  felm(I(100*lodage4)~year|country|0|year,
                     data=seg_age4[(seg_age4$country %in% "France" & seg_age4$year==1993)==F
                                   & (seg_age4$country %in% "Denmark" & seg_age4$year==1994)==F
                                   ,])

m_oocc1_all <-  felm(I(100*lodocc1)~year|country|0|year,
                     data=seg_occ1[(seg_occ1$country %in% "France" & seg_occ1$year==1993)==F
                                   & (seg_occ1$country %in% "Denmark" & seg_occ1$year==1994)==F
                                   ,])

m_oocc3_all <-  felm(I(100*lodocc3)~year|country|0|year,
                     data=seg_occ3[(seg_occ3$country %in% "France" & seg_occ3$year==1993)==F
                                   & (seg_occ3$country %in% "Denmark" & seg_occ3$year==1994)==F
                                   ,])



#table S5.1B
screenreg(list(m_eo55,m_eo4545,m_mig_eo,m_gdr_eo,m_oage1,m_oage4,m_oocc1,m_oocc3),
          custom.model.names = c("m_eo55","m_eo4545","m_mig_eo","m_gdr_eo","m_oage1","m_oage4","m_oocc1","m_oocc3"),
          digits=1,stars=c(0.1,0.05,0.01))
htmlreg(list(m_eo55,m_eo4545,m_mig_eo,m_gdr_eo,m_oage1,m_oage4,m_oocc1,m_oocc3),
        custom.model.names = c("m_eo55","m_eo4545","m_mig_eo","m_gdr_eo","m_oage1","m_oage4","m_oocc1","m_oocc3"),
        digits=1,stars=c(0.1,0.05,0.01),file="Tables/Table S5.1B.html")

#table S5.1C
screenreg(list(m_eo55_all,m_eo4545_all,m_mig_eo_all,m_gdr_eo_all,m_oage1_all,m_oage4_all,m_oocc1_all,m_oocc3_all),
          custom.model.names = c("m_eo55_all","m_eo4545_all","m_mig_eo_all","m_gdr_eo_all","m_oage1_all","m_oage4_all","m_oocc1_all","m_oocc3_all"),
          digits=1,stars=c(0.1,0.05,0.01))
htmlreg(list(m_eo55_all,m_eo4545_all,m_mig_eo_all,m_gdr_eo_all,m_oage1_all,m_oage4_all,m_oocc1_all,m_oocc3_all),
        custom.model.names = c("m_eo55_all","m_eo4545_all","m_mig_eo_all","m_gdr_eo_all","m_oage1_all","m_oage4_all","m_oocc1_all","m_oocc3_all"),
        digits=1,stars=c(0.1,0.05,0.01),file="Tables/Table S5.1C.html")





